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Preface 


This  report  describes  the  results  of  a  three-dimensional  Finite  Element 
seepage  analysis  of  Cerrillos  Dam  conducted  during  the  period  November 
1992  to  July  1993.  This  dam,  located  near  Ponce,  Puerto  Rico,  was  recently 
constructed  and  the  observed  foundation  seepage  was  in  excess  of  the  design 
estimate.  This  study  was  undertaken  as  an  effort  to  evaluate  the  reasonable¬ 
ness  of  the  observed  seepage  discharge  rates.  The  study  was  requested  and 
funded  by  the  U.S.  Army  Engineer  District,  Jacksonville  (CESAJ-EN-G). 
Technical  Monitor  was  Mr.  Paul  Shafer,  Chief,  Geotechnical  Branch. 

Dr.  Edward  E.  Middleton  was  Chief,  Engineering  Division,  and  COL 
Terrence  Salt  was  Commander  of  the  Jacksonville  District. 

The  study  was  conducted  and  the  report  was  prepared  by  Mr.  John  B. 
Palmerton,  Rock  Mechanics  Branch,  Soil  and  Rock  Mechanics  Division, 
Geotechnical  Laboratory,  U.S.  Army  Engineer  Waterways  Experiment  Station 
(CEWES-GS-R).  The  study  was  conducted  under  the  direct  supervision  of 
Mr.  Jerry  S.  Huie,  Chief,  Rock  Mechanics  Branch,  and  under  the  general 
supervision  of  Dr.  Don  C.  Banks,  Chief,  Soil  and  Rock  Mechanics  Division, 
and  Dr.  William  F.  Marcuson  m.  Director,  Geotechnical  Laboratory. 

Ms  Kristine  McNamara,  Soil  Testing  Facility  (CEWES-GS-GD)  prepared  the 
graphical  presentation  portions  of  this  report. 

At  the  time  of  publication  of  this  report,  Director  of  WES  was 
Dr.  Robert  W.  Whalin.  Commander  was  COL  Bruce  K.  Howard,  EN. 


Conversion  Factors, 
Non-SI  to  SI  Units  of 
Measurement 


1  Introduction 


Cerrillos  Dam  is  located  near  Ponce,  Puerto  Rico,  on  the  Cerrillos  River. 
The  dam  is  approximately  323  ft  high  with  a  crest  length  of  1,555  ft.  Con¬ 
struction  of  this  dam  was  completed  in  1992.  Figure  1  is  a  photo  taken  in 
March  1992  which  shows  the  completed  structure  as  seen  looking  east  from 
the  top  of  the  right  abutment.  The  upstream  curvature  of  the  dam  is  readily 
apparent  in  the  photo  as  is  the  inclined  outlet  structure.  The  outlet  structure  is 
not  integral  to  the  embankment,  but  is  founded  on  the  sloping  face  of  a  small 
ridge  which  extends  toward  the  pool.  The  left  valley  wall  in  the  vicinity  of 
the  outlet  structure  was  scaled  and  cleared  following  a  landslide  which 
occurred  during  the  foundation  excavation  of  the  outlet  structure.  The  dam 
consists  of  a  central  clay  core,  a  grout  curtain  extending  to  a  depth  of  200  ft, 
and  upstream  and  downstream  rockfill  shells  with  appropriate  transition  zones. 
The  first  filling  of  the  dam  occurred  in  late  1992  and  considerable  seepage  (as 
intercepted  by  a  toe  drain  system  in  the  downstream  shell)  was  observed. 
Analysis  of  the  toe  drain  observations  indicates  that  approximately  4  cubic  ft 
per  second  (cfs)  was  being  collected  by  the  drain  system  at  a  pool  depth  of 
approximately  200  ft  (el  495  ft).  The  planned  maximum  conservation  pool 
will  be  281  ft  deep  (el  573  ft).  The  design  seepage  estimate  was  1.0  cfs  at  the 
maximum  conservation  pool. 

The  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES)  was 
requested  by  the  U.S.  Army  Corps  of  Engineer  District,  Jacksonville,  to 
undertake  a  three-dimensional  (3-D)  seepage  analysis  of  Cerrillos  Dam  in 
November  1992.  This  report  documents  the  results  of  the  analyses. 


Purpose 

The  purpose  of  this  study  is  to  conduct  a  thorough  analysis  of  the  seepage 
regime  at  the  Project  site  and  to  provide  a  tool  to  evaluate  the  effectiveness  of 
different  concepts  for  seepage  control  at  the  site.  It  was  specifically  requested 
by  the  District  that  a  3-D  finite  element  (FE)  analysis  be  conducted. 
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Background 


In  November  1992,  the  author  visited  the  Jacksonville  District  offices  and 
was  briefed  by  District  staff  engineers  and  geologists  as  to  the  pertinent  details 
relating  to  the  observed  seepage.  At  the  time  of  this  briefing,  field  evidence 
indicated  that  the  majority  of  the  observed  seepage  was  coming  from  the 
direction  of  the  left  abutment.  A  strong  feeling  existed  among  the  District 
staff  that  the  seepage  was  more  than  likely  due  to  the  presence  of  a  geologic 
structural  unit  known  as  the  Ridge  (Forming)  Limestone.  Extensive  geologic 
studies  had  previously  been  conducted  and  their  results  incorporated  in  the 
various  Project  Design  Memoranda,  particularly  the  Feature  Design 
Memorandum  No.  171  (DM-17)  which  was  issued  in  1983.  DM-17  presents 
(among  other  items)  detailed  geologic  structure  maps,  the  seepage  analysis 
calculations  for  the  designed  project  (2-D),  and  results  from  field  tests  to 
determine  the  permeability  of  the  foundation  materials.  It  was  decided  during 
the  November  meeting  to  conduct  a  3-D  FE  seepage  analysis  which  specifi¬ 
cally  included  the  effects  of  the  Ridge  Limestone  unit. 


*  In  this  report,  elevations  are  in  feet  National  Geodetic  Vertical  Datum  (NGVD). 
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Pertinent  Features  at  Cerriilos  Dam  Influencing 
Seepage 


The  geologic  structure  in  the  vicinity  of  the  dam  is  characterized  by  steeply 
dipping  planar  and  parallel  units  of  limestone,  siltstone,  and  tuff  which  under¬ 
lay  a  very  rugged  terrain.  Surficial  soils  are  either  non-existent  or  are  very 
thin.  The  bedding  strikes  N  50-60°  E  and  dips  40-50°  SE.  Permeability 
coefficients  within  these  units  show  a  wide  variation  with  the  Ridge  Limestone 
exhibiting  the  highest  (by  far)  permeability  (see  page  5-3  of  DM-17).  The 
Ridge  Limestone  unit  outcrops  on  the  left  valley  wall  in  the  vicinity  of  the 
location  of  the  dam’s  clay  core  and  the  left  abutment  interface  as  shown  in 
Figure  2  (taken  from  Plate  10-22  of  DM-17).  At  the  base  of  the  valley,  the 
Ridge  Limestone  unit  is  exposed  beneath  the  downstream  rockfill  shell  of  the 
dam  as  shown  in  Figure  3  (taken  from  Plate  10-19  of  DM-17).  During  con¬ 
struction,  a  grout  curtain  to  reduce  seepage  was  placed  by  injecting  grout  in 
drill-holes  to  a  depth  of  200  ft  under  the  full  length  of  the  impervious  clay 
core  contact  with  the  valley  floor  and  abutments.  The  pertinent  data  for  esti¬ 
mating  the  seepage  discharge  quantities  during  the  design  phase  are  shown  in 
Figure  4.  This  analysis  was  accomplished  by  constructing  a  flow  net  (using 
transformed  sections  to  acount  for  the  variation  of  permeability  with  depth). 
This  analysis  did  not  specifically  include  the  Ridge  Limestone  unit  nor  did  it 
include  the  seepage  through  the  abutments  above  the  valley  floor.  This  type 
of  flow  net  analysis  is  the  typical  procedure  for  Corps  of  Engineers  dams.  In 
most  c  ises  the  seepage  around  dam  abutments  is  conservatively  (since  the 
pressure  head  decreases  in  the  abutment  zones)  accounted  for  by  assuming  a 
crest  length  which  is  in  excess  of  the  width  of  the  valley  floor.  However,  if 
the  permeability  of  the  materials  within  the  abutments  (i.e.,  the  material  above 
the  valley  floor)  differs  significantly  from  those  beneath  the  valley  floor,  a 
different  type  analysis  would  be  indicated.  For  Cerriilos  Dam,  the  failure  to 
specifically  account  for  the  presence  of  the  valley  crossing  Ridge  Limestone 
unit  apparently  led  to  a  low  estimate  for  the  quantity  of  seepage.  The  topog¬ 
raphy  of  the  region  in  the  vicinity  of  the  dam  is  shown  in  Figure  5.  The 
position  (outcrop)  of  the  Ridge  Limestone  unit  and  the  location  of  the  dam, 
spillway,  and  intake  tower  is  shown  in  Figure  6. 

The  3-D  FE  method  analyses  required  some  simplifications  of  the  actual 
geometry.  In  particular,  the  valley  floor  and  walls  were  assumed  to  be 
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uniform  and  straight  in  the  vicinity  of  the  dam  and  the  cross-valley  curvature 
of  the  dam  was  eliminated.  These  assumptions  are  illustrated  in  Figure  7 
which  depicts  the  relationship  of  the  dam,  grout  curtain,  and  the  Ridge  Lime¬ 
stone.  Figure  8,  in  which  the  rockfill  shell  is  depicted  as  being  transparent, 
shows  the  exposure  of  the  Ridge  Limestone  upstream  of  the  dam  and  beneath 
the  downstream  rockfill  shell.  This  figure  also  shows  the  location  of  the  tuff 
and  siltstone  units  which  parallel  the  Ridge  Limestone.  The  location  of  the 
seepage  collection  system  is  shown  in  Figure  9.  A  cut-away  view  showing 
various  zones  of  the  dam  and  foundation  is  shown  in  Figure  10.  A  "birds- 
eye"  view  is  shown  in  Figure  11;  again,  the  shells  of  the  dam  are  depicted  as 
being  transparent.  A  cross  section  taken  parallel  to  the  valley  and  midway  up 
the  left  valley  wall  (looking  downstream)  is  shown  in  Figure  12.  The  Ridge 
Limestone  unit  is  positioned  such  that  water  may  percolate  from  the  pool  into 
that  rock  unit  along  the  left  valley  wall  (just  upstream  of  the  dam  and  within 
the  dam’s  permeable  upstream  rockfill  shell)  and  be  carried  in  a  relatively 
unimpeded  fashion  (compared  with  the  surrounding  rock)  under  the  grout 
curtain  (which  also  intercepts  a  portion  of  the  Ridge  Limestone)  and  then 
emerges  from  beneath  the  downstream  rockfill  shell  and  from  the  lower  eleva¬ 
tions  of  the  downstream  left  abutment  face.  Seepage  water  which  is  tran¬ 
sported  in  this  manner  would  be  collected  by  the  toe  drain  system.  This  mode 
by  which  seepage  water  is  transported  is  also  indicated  in  Figure  13  (in  which 
the  geologic  materials  downstream  of  the  Ridge  Limestone  are  removed  from 
the  drawing). 


Finite  Element  Mesh  Generation 

A  very  recently  developed  3-D  FE  method  seepage  computer  program 
which  is  operational  on  the  Cray  Y-MP  supercomputer  at  WES  was  used  for 
die  analyses  of  Cerrillos  Dam.  A  complete  explanation  of  the  theory  and 
application  of  this  computer  program  is  given  by  Tracy  (1991).  To  use  the 
program,  a  3-D  mesh  (consisting  of  node  points  and  elements)  must  be  pre¬ 
pared.  Each  element  is  assigned  a  material  type  which  is  characterized  by  a 
permeability  tensor.  That  is,  different  Darcy  law  coefficients  of  permeability 
may  be  prescribed  in  each  of  the  three  coordinate  directions.  Boundary  piezo¬ 
metric  heads  and  boundary  flow  quantities  or  flow  velocities  can  be  specified 
as  part  of  the  input  data.  This  FE  formulation  can  simulate  confined  or 
unconfined  flow  regimes.  For  unconfined  flow  problems  (such  as  Cerrillos 
Dam)  the  location  of  the  phreatic  surface  is  established  by  an  iterative  process. 
The  output  of  the  program  consists  of  the  total  seepage  discharge  quantity,  the 
piezometric  head  at  each  node  point  and  the  flow  velocities.  Both  the  FE 
method  and  the  graphical  flow  net  technique  seek  solutions  to  identical  mathe¬ 
matical  formulations  for  seepage  flow  (i.e.,  Laplace’s  equation  and  Darcy’s 
law);  therefore,  either  method  will  yield,  excluding  computational  errors, 
identical  results.  The  FE  method  may  be  used  to  simulate  a  much  larger 
range  of  seepage  applications  since  flow  net  solutions  are  limited  to  2-D  appli¬ 
cations  which  must  be  reduced  (by  scaling  transformations)  to  equivalent 
isotropic  (uniform  permeability  in  all  direction)  situations. 
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To  perform  the  many  analyses  conducted  for  Cerrillos  Dam  during  this 
study,  the  preparation  of  a  3-D  FE  mesh  was  the  first  step.  Since  3-D  FE 
analysis  requires  so  many  nodes  and  elements  to  describe  even  a  very  simple 
seepage  problem,  it  was  necessary  to  write  a  mesh  generator  computer  pro¬ 
gram  (called  MESH).  The  MESH  program  (a  listing  is  given  in  Appendix  A) 
was  specifically  tailored  to  account  for  the  following  pertinent  details  of 
Cerrillos  Dam: 

a.  The  geologic  units,  especially  the  Ridge  Limestone  unit  must  strike  and 
dip  as  discussed  in  earlier. 

b.  The  elements  located  in  the  vertical  dam  clay  core  must  be  identified. 

c.  The  elements  located  within  the  grout  curtain  must  be  identified. 

d.  The  position  of  the  upstream  and  downstream  dam  shells  must  be  identi¬ 
fied.  Because  of  their  very  high  permeability,  the  dam  rockfill  shells 
were  not  included  in  the  FE  mesh.  However,  the  position  of  the  shells 
were  needed  to  simulate  the  effect  of  corrective  treatment  (i.e. 
shotcreting). 

e.  Depending  on  the  height  of  pool  being  considered,  the  node  points  for 
which  boundary  piezometric  heads  are  to  be  applied  must  be  identified. 

The  complete  details  of  the  automatic  mesh  generator  are  too  lengthy  for 
inclusion  herein.  In  brief,  a  box-like  (see  Figure  14a)  mesh  with  a  pair  of 
each  element’s  faces  slanting  parallel  to  the  strike  and  dip  of  the  geologic  units 
(specified  to  be  striking  at  N  55°  E  and  dipping  50.7°  SE)  was  first  gener¬ 
ated.  The  height  of  the  box-like  grid  is  dependent  upon  the  specified  pool 
height.  The  North  axis  was  taken  to  be  in  the  positive  y  direction,  the  East 
axis  is  in  the  positive  x  direction,  and  z  is  positive  upward.  The  above  values 
for  the  strike  and  dip  leads  to  an  apparent  dip  of  35°  E  with  the  x-z  (or  cross 
valley)  plane  (which  is  also  the  same  angle  that  the  valley  walls  make  with  the 
horizontal  plane),  and  an  apparent  dip  of  45°  S  with  the  y-z  plane  (or  longitu¬ 
dinal  plane).  Following  the  generation  of  the  box-like  mesh,  the  valley  is 
created  by  locating  the  elements  whose  centers  lie  aoove  the  valley  floor  and 
between  die  two  valley  walls  (see  Figure  14b).  All  of  those  elements,  except¬ 
ing  those  within  the  dam’s  clay  core,  are  then  removed.  The  nodes  and  ele¬ 
ments  are  renumbered  so  that  only  those  nodes  and  elements  which  are  still 
present  (i.e.,  not  removed)  shall  be  used  in  the  FE  analysis.  Those  elements 
which  lie  between  the  planes  delineating  the  grout  curtain  and  the  dam  core 
are  then  associated  with  the  material  type  for  those  zones.  Finally,  depending 
on  the  specified  (for  a  particular  problem)  height  of  pool,  the  nodes  which  lie 
on  the  valley  floor  and  valley  walls  upstream  of  the  dam  and  on  the  upstream 
face  of  the  clay  core  are  identified  for  specification  of  the  boundary  piezomet¬ 
ric  head  (e.g.,  200  ft).  In  a  similar  fashion,  the  boundary  piezometric  head 
on  the  nodes  on  the  downstream  valley  and  downstream  face  of  the  clay  core 
are  set  to  the  tail  water  head  (zero). 
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The  overall  dimensions  of  the  box-like  mesh  is  3,141.6  ft  in  the  x,  or  cross 
valley  dimension,  2,000  ft  in  the  y,  or  longitudinal  direction,  and  a  maximum 
height,  z  of  1,200  ft  as  shown  in  Figure  15.  (Fcr  each  case  for  which  an  FE 
analysis  was  made,  the  height  of  the  box  was  set  at  the  variable  elevation  of 
the  pool).  The  width  of  the  valley  floor  (which  is  assumed  to  be  a  horizontal 
plane)  is  428  ft.  The  valley  walls  rise  to  a  maximum  height  of  350  ft  at  an 
angle  of  35°  to  the  horizontal.  For  the  case  of  a  350-ft  pool  height,  the  num¬ 
ber  of  (final)  generated  elements  is  8,283,  with  10,810  nodes;  for  a  100-ft 
high  pool,  7,340  elements  and  9,405  nodes  are  generated. 


Material  Properties  and  Model  Calibration 

As  previously  stated  in  paragraph  5,  a  good  amount  of  field  permeability 
testing  was  conducted  in  years  past  for  the  Cerrillos  Dam  project.  These  field 
tests  led  to  selection  of  the  permeability  coefficients  listed  in  Figure  4.  The  k 
(Darcy  permeability  coefficients)  values  ranged  from  1.4  x  10'3  cm/sec  near 
the  ground  surface  to  1.4  x  10"4  cm/sec  at  depths  below  the  surface  greater 
than  200  ft.  Water  pressure  injection  tests  conducted  at  the  site  indicated 
permeability  values  ranging  from  5.1  x  10'3  cm/sec  in  the  Ridge  Limestone  to 
2.5  x  10'5  cm/sec  at  greater  depths.  A  decrease  in  permeability  with  depth  is 
almost  always  observed  at  most  sites.  This  decrease  is  due  to  overburden 
caused  stresses  causing  cracks  and  rock  joints  to  be  more  tightly  closed  at 
depth.  The  2-D  case  (with  the  k  values  as  listed)  shown  in  Figure  4  was 
simulated  with  the  3-D  FE  program  by  setting  the  permeability  coefficients 
beyond  the  base  width  of  the  valley  floor  to  a  very  low  value  (1  x  10'12 
cm/sec)  and  then  multiplying  the  computed  seepage  discharge  by  the  ratio  of 
the  length  of  the  hand-calculated  section  (1,555  ft)  to  the  width  of  the  valley 
floor  in  the  3-D  model  (428  ft).  As  desired,  a  result  very  near  to  the  seepage 
discharge  of  1  cfs  noted  in  Figure  4  was  computed  by  the  FE  program.  Fur¬ 
ther  2-D  simulations  were  made  to  determine  which  uniform  rock  permeability 
value  would  also  yield  a  seepage  discharge  of  1  cfs.  A  k  of  1.22  x  10"4 
cm/sec  assigned  to  all  rock  below  the  valley  floor  yielded  1  cfs.  Because  of 
the  steeply  dipping  beds  and  the  often  apparent  decrease  in  permeability  with 
depth,  no  attempt  was  made  to  consider  anisotropic  permeability  properties. 

A  rational  method  for  assigning  principle  permeability  tensor  directions  is  too 
difficult. 


3-D  Finite  Element  Computations 

Even  though  the  3-D  FE  program  (employed  for  2-D  simulations)  yielded  a 
close  comparison  to  the  hand  calculation,  the  issue  of  what  k  values  to  assign 
to  the  rock  above  the  valley  floor  was  not  yet  resolved.  Various  trial  com¬ 
puter  simulations  for  which  various  combinations  of  k  values  for  the  Ridge 
Limestone,  the  surrounding  rock  formation,  and  the  permeability  of  the  grout 
curtain  were  conducted.  After  an  examination  and  discussion  of  these  simula¬ 
tions  with  District  engineers  and  geologists,  the  set  of  permeability  coefficients 
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listed  in  Table  1  was  selected.  The  delineation  of  the  zones  is  shown  in 
Figure  15. 

The  combination  of  permeability  coefficients  listed  in  Table  1  yielded  an 
excellent  estimate  for  the  observed  seepage  discharges  at  the  project  for  a  pool 
height  of  about  200  ft.  The  predicted  seepage  discharge  versus  pool  elevation 
is  shown  in  Figure  16  0ower  <r.  v^).  If  the  pool  elevation  is  permitted  to 
continue  its  rise  (without  furth* »  remedial  efforts  to  control  the  seepage),  a 
seepage  discharge  of  approximately  6.5  cfs  (at  the  maximum  conservation 
pool)  is  predicted.  Data  which  are  gathered  as  the  pool  rises  in  the  future 
may  also  be  used  to  modify  the  material  properties  in  order  to  refine  the 
predictions.  Table  1  indicates  that  the  rock  below  el  105  (400  ft  below  the 
valley  floor)  was  considered  to  be  impermeable.  The  upper  curve  in  Fig¬ 
ure  16  shows  the  effect  of  assigning  a  k  value  of  1  x  1(P  cm/sec  to  the  rock 
below  el  105  ft.  Only  a  small  increase  in  seepage  discharge  is  predicted. 
Impermeable  lower  rock  zones  were  desirable  to  maintain  compatibility  with 
the  2-D  flow  net  calculations.  Typically,  each  iterative  3-D  FE  analysis 
required  a  total  of  2  to  3  min  of  processor  time  (and  10  to  15  min  of  wall 
clock  time)  for  completion  on  the  supercomputer.  The  MESH  grid  generator 
program  required  approximately  12  sec  of  supercomputer  processor  time  to 
generate  about  8,000  elements  and  11,000  nodes. 


Table  1 

Grout  curtain  (200-ft  depth  below  clay  core) 

1.0  x  10  ®  cm/sec 

Impermeable  clay  core  of  dam 

6.4  x  10'7  cm/sec 

Ridge  Limestone  unit 

5.0  x  10‘3  cm/sec 

All  other  foundation  rock  units  above  el  -105  ft1 

1 .0  x  1 0"4  cm/sec 

All  foundation  rock  units  below  el  -105  ft 

Impermeable 

1  El  105  is  400  ft  below  the  valley  floor  which  is  at  el  295. 

The  3-D  FE  analyses  were  not  conducted  for  the  purpose  of  predicting 
seepage  discharges  at  this  specific  site.  Prediction  is  not  the  issue  since  the 
answer  is  already  known.  Further,  since  the  seepage  discharges  are  generally 
linear  with  respect  to  pool  height,  the  simple  extrapolation  of  field  observa¬ 
tions  would  provide  as  good  (or  better)  a  predictive  capability  as  the  computer 
model.  Instead,  the  usefulness  of  this  study  has  been  the  demonstration  that 
the  observed  higher  than  estimated  seepage  discharges  at  Cerrillos  Dam  have 
been  replicated  in  a  rational  manner  that  did  not  require  any  major  modifi¬ 
cation  of  the  understanding  of  the  rock  permeability  properties  which  existed 
during  the  design  phase.  It  is  doubtful  that  2-D  (whether  FE  of  flow  net) 
analysis  could  have  rationally  predicted  the  higher  seepage  discharges.  The 
complicated  geologic  conditions  necessitate  a  3-D  analysis.  Had  a  3-D  FE 
seepage  analysis  of  Cerrillos  Dam  been  attempted  during  the  design  phases  of 
the  project,  it  is  questionable  if  the  now-measured  seepage  discharges  would 
have  been  exactly  predicted.  However,  since  the  Ridge  Limestone  had  been 
identified  as  a  very  permeable  geologic  unit,  any  reasonable  3-D  analysis  in 
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which  different  permeability  combinations  were  varied  would  certainly  have 
served  to  alert  the  designers  of  a  potential  problem.  The  foundation  perme¬ 
ability  model  presented  in  Table  1  presumes  that  the  Ridge  Limestone  is 
SO  times  more  permeable  that  the  other  foundation  rock.  Field  permeability 
tests  bear  out  the  reasonableness  of  the  modeled  permeability  coefficient 
values. 

A  useful  and  important  application  of  3-D  FE  modeling  is  the  evaluation  of 
the  effectiveness  of  engineered  modifications  for  seepage  control.  For  exam¬ 
ple,  it  has  been  suggested  that  the  placement  of  a  layer  of  (impermeable)  shot- 
crete  on  the  upstream  valley  wall  exposures  of  the  Ridge  Limestone  might  be 
effective  in  lowering  the  seepage  discharges.  However,  during  the  construc¬ 
tion  of  the  dam,  the  upstream  rockfUl  shell  was  placed  against  the  exposed 
Ridge  Limestone  on  the  left  abutment.  Therefore,  it  will  not  be  possible  to 
shotcrete  the  Ridge  Limestone  buried  beneath  the  upstream  shell  without  the 
shell’s  removal.  It  was  fairly  simple  to  modify  the  MESH  program  so  that 
the  layer  of  elements  within  the  Ridge  Limestone  adjacent  to  the  pool  could  be 
made  more  impermeable  (i.e.  changed  to  k=6.4  x  10'7  cm/sec).  The  simula¬ 
tions  of  shotcreting  the  Ridge  Limestone  (a)  only  above  the  upstream  shell, 

(b)  only  below  the  upstream  shell,  and  (c)  the  simulation  of  shotcreting  the 
entire  left  valley  wall  upstream  exposure  were  undertaken.  Figure  17  shows 
the  result  of  these  simulations.  The  upper  curve  shows  the  predicted  seepage 
discharge  for  no  placement  of  shotcrete  (this  is  also  the  lower  curve  shown  in 
Figure  16).  The  next  lower  curve  shows  the  seepage  discharge  prediction  for 
placement  of  shotcrete  only  on  the  now  visible  (and  accessible)  upstream 
Ridge  Limestone  face  above  the  upstream  shell.  The  third  lower  curve  shows 
the  seepage  discharge  for  placement  of  shot^ete  (or  perhaps  grout)  only  on 
the  exposure  of  Ridge  Limestone  now  buried  by  the  upstream  rockfUl  shell. 
The  fourth  lower  curve  is  for  shotcrete  placement  on  all  of  the  upstream 
exposed  valley  wall  limestone.  It  is  apparent  from  these  curves  that  partial 
shotcreting  will  not  be  nearly  as  effective  as  complete  coverage  of  the  highly 
permeable  Ridge  Limestone  unit.  As  the  pool  is  raised,  the  portions  of  the 
untreated  limestone  beneath  the  shell  will  still  cause  considerable  leakage 
because  the  water  pressure  head  along  those  untreated  faces  will  continue  to 
rise.  The  lower  curve  in  Figure  17  is  indicative  of  the  best  that  could  be 
expected  with  any  treatment  of  the  Ridge  Limestone.  This  curve  was  obtained 
by  presuming  that  the  permeability  of  the  Ridge  Limestone  is  the  same  as  the 
surrounding  rock  (i.e.  k=  1  x  10~*  cm/sec).  Also  notice  that  the  predicted 
seepage  discharge  for  this  "best  case"  scenario  is  approximately  1.0  cfs  (the 
2-D  flow  net  "design"  estimate)  at  the  maximum  conservation  pool  elevation. 

The  effect  of  deepening  the  grout  curtain  to  lessen  the  amount  of  seepage 
discharge  is  shown  in  Figure  18.  It  is  seen  from  this  plot  that  an  increase  in 
depth  of  100  ft  (to  a  total  depth  of  300  ft)  would  reduce  the  seepage  quantity 
to  S.3  cfs  (from  the  estimated  seepage  quantity  of  6.5  cfs  at  the  conservation 
pool  elevation),  an  18  percent  improvement.  Doubling  the  depth  of  the  grout 
curtain  to  400  ft  would  reduce  the  seepage  quantity  to  4.5  cfs,  a  31  percent 
improvement.  The  flattening  of  the  curve  at  larger  grout  curtain  depths  is  not 
favorable  to  this  method  of  seepage  reduction. 
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The  3-D  FE  program  produces  much  more  output  information  than  simply 
the  total  seepage  discharge.  Other  types  of  output  include  the  piezometric 
head  at  each  node  point,  the  position  of  the  phreatic  surface,  and  the  flow 
velocities  within  each  element.  For  example,  a  3-D  plot  of  the  piezometric 
heads  is  shown  in  Figure  19.  A  cross  section  along  the  valley  center  shown 
in  Figure  20  illustrates  the  high  head  gradient  across  the  dam  core  and  grout 
curtain.  Flow  velocity  plots  are  shown  in  Figures  21,  22  and  23.  These  plots 
indicate  the  absolute  magnitude  of  the  velocity  without  regard  to  direction  of 
flow.  Figure  21  shows  that  the  region  of  higher  flow  velocities  are  concen¬ 
trated  in  the  zones  which  comprise  the  Ridge  Limestone.  Figure  22  shows  the 
concentration  of  seepage  velocity  along  a  vertical  cross  section  just  down¬ 
stream  of  the  dam  core;  again,  the  higher  velocities  are  concentrated  within 
the  Ridge  Limestone.  The  seepage  velocities  along  a  horizontal  section  taken 
immediately  above  the  base  of  the  dam  are  shown  in  Figure  23.  Because  of 
the  large  number  of  nodes  and  elements  required  for  practical  problems,  it 
was  necessary  to  write  additional  computer  programs  to  process  the  output 
data  file  generated  by  the  3-D  FE  program.  Four  computer  programs,  named 
HEADS,  PREFLOW,  FLOWS,  and  FLOVEC  (see  Appendix  A)  were  written 
during  the  course  of  this  study  to  process  the  FE  output  data  files.  These 
programs  were  written  in  FORTRAN  for  a  personal  computer.  The  FE  out¬ 
put  file(s)  (about  3  Mbytes  in  size)  were  retrieved  from  the  supercomputer 
over  the  WES  fiber  optic  network  and  stored  on  the  personal  computer.  The 
HEADS  program  examined  the  FE  output  file(s)  and  sorted  and  stored  the 
data  into  smaller  sub-files.  Each  sub-file  contained  the  piezometric  head  val¬ 
ues  on  numerous  planar  sections  passing  through  the  FE  mesh.  These 
sub-files  were  then  processed  by  a  contouring  program  named  McCON 
(Palmerton  1992).  A  similar  function  was  performed  by  the  PREFLOW  and 
FLOWS  programs  for  the  output  of  seepage  velocities  at  the  element  centers. 
The  FLOVEC  program  prepared  graphical  plots  of  the  flow  velocity  direction 
vectors. 

Examples  of  the  contour  plots  are  shown  in  Figures  24  through  26.  The 
contours  of  piezometric  head  on  a  horizontal  plane  at  el  45  ft  (250  ft  below 
the  base  of  the  dam)  caused  by  a  pool  depth  of  200  ft  (el  495  ft)  are  shown  in 
Figure  24.  The  contour  values  are  in  terms  of  elevation.  Notice  the  high 
head  gradient  near  the  center  of  the  plot  where  the  water  flow  is  restricted  by 
the  grout  curtain.  Also  notice  that  the  heads  are  somewhat  higher  (on  the 
order  of  15-20  ft)  toward  the  left  abutment  (right  side  of  plot).  The  higher 
heads  are  due  to  the  presence  of  the  more  permeable  Ridge  Limestone  in  the 
left  valley  wall.  The  piezometric  head  contours  at  el  245  (50  ft  below  the 
base  of  the  dam)  are  shown  in  Figure  25.  Piezometric  head  contours  along  a 
vertical  longitudinal  plane  located  at  x  =  1,713  ft  (see  Figure  6  for  orienta¬ 
tion)  are  shown  in  Figure  26.  The  piezometric  heads  decrease  quite  rapidly 
across  the  clay  core  and  grout  curtain. 

Vector  plots  of  the  computed  direction  of  flow  along  various  sections  are 
shown  in  Figures  27  through  35.  Plots  of  this  type  are  valuable  aids  in 
visualizing  the  character  of  the  seepage  paths.  Figures  27,  28,  and  29  show 
the  direction  of  the  flow  on  horizontal  planes  situated  at  el  70  (225  ft  below 
the  base  of  the  dam  and  25  ft  below  the  grout  curtain),  el  270  (25  ft  below  the 
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base  of  the  dam),  and  el  370  (75  ft  above  the  base  of  the  dam).  The  vectors 
from  the  larger  circles  denote  the  direction  of  the  flow  in  the  Ridge  Lime¬ 
stone.  As  the  elevation  of  the  sections  increase,  the  position  of  the  Ridge 
Limestone  (due  to  the  dip  angle  of  the  unit)  apparently  migrates  toward  the 
upper  left  on  the  plots.  The  small  circles  and  large  solid  circles  indicate  that 
the  component  of  the  velocity  is  toward  the  eye;  the  Vs"  indicate  the  oppo¬ 
site.  All  of  these  figures  give  the  result  for  a  specified  pool  height  of  200  ft 
(and  a  total  seepage  discharge  of  4.3  cfs  -  the  case  for  no  shotcrete  place¬ 
ment).  The  data  points,  which  appear  to  be  missing  in  Figure  29,  are  the 
result  of  the  3-D  FE  program  collapsing  the  mesh  so  that  no  points  lie  above 
the  phreatic  surface.  The  phreatic  surface  is  also  somewhat  higher  in  eleva¬ 
tion  on  the  left  side  of  the  valley  (right  side  of  plot)  as  will  be  discussed 
below. 

Vector  flow  plots  on  vertical  longitudinal  sections  located  at  x=  1,606, 
2,034,  and  2,320  ft  (see  Figure  6  for  orientation)  are  shown  in  Figures  30, 

31,  and  32.  These  sections  are  taken  near  the  center  of  the  valley  floor, 
mid-way  up  the  left  valley  wall  at  approximately  the  pool  elevation,  and  paral¬ 
lel  to  the  left  ridge  some  35  ft  beyond  the  crest.  Again,  the  large  circles  indi¬ 
cate  the  location  of  the  Ridge  Limestone  on  each  of  the  sections.  The  200-ft 
pool  height  is  indicated  on  the  plots.  It  is  interesting  to  note  that  Figure  30 
shows  the  direction  of  flow  within  the  Ridge  Limestone  to  be  upstream  as  it 
emerges  upward  from  the  valley  floor.  Again,  the  small  circles  and  the  large 
solid  circles  indicate  that  the  (cross  valley)  component  of  the  flow  is  toward 
the  eye,  and  the  Vs"  indicate  that  the  flow  is  opposite,  or  away  from  the 
eye. 

Vector  flow  plots  on  vertical  cross  valley  sections  are  shown  in  Figures  33, 
34,  and  35.  The  section  shown  in  Figure  33  is  taken  approximately  250  ft 
downstream  of  the  axis  of  the  dam.  Figure  34  is  taken  parallel  and  through 
the  axis  of  the  dam  and  Figure  35  is  taken  250  ft  upstream.  Downstream  of 
the  dam  (Figure  33),  the  flows  are  directed  toward  the  valley;  upstream  (Fig¬ 
ure  34),  the  flows  are  generally  directed  away  from  the  valley  with  some 
exceptions  in  the  Ridge  Limestone  (large  circles).  For  die  section  along  the 
dam  axis  (Figure  35),  the  flows  are  still  directed  outward  from  the  valley  on 
the  right  valley  side  (left  side  of  the  plot),  but  trading  to  flow  toward  the 
valley  on  the  left  valley  wall.  These  plots  all  indicate  die  pool  elevation  and 
the  location  of  the  grout  curtain. 

The  3-D  FE  computer  program  computes  the  location  of  the  phreatic  sur¬ 
face  and  adjusts  the  mesh  so  that  the  portions  lying  above  this  surface  collapse 
to  the  surface.  For  an  unconfined  seepage  analysis,  the  phreatic  surface  is 
defined  to  be  that  surface  on  which  the  elevation  head  is  equal  to  the  piezo¬ 
metric  head.  The  output  data  from  the  3-D  FE  analysis  (for  a  pool  height  of 
200  ft)  were  examined  by  a  program  entitled  PHREAT  (see  Appendix  A) 
written  during  this  study  and  the  points  (and  their  coordinate  location)  which 
were  identified  as  being  on  the  phreatic  surface  were  extracted.  These  data 
points  were  then  processed  by  the  McCON  contouring  program  and  contours 
of  the  phreatic  surface  were  plotted  as  shown  in  Figure  36.  The  McCON 
contouring  program  may  also  be  used  to  prepare  profiles  on  multiple  sections 
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perpendicular  to  the  plane  of  the  contour  plots.  Figure  37  shows  the  profiles 
of  two  sections  taken  downstream  of  the  dam  axis  and  two  sections  taken 
upstream  (see  Figure  21  to  locate  the  section  numbers).  Notice  in  Figure  37 
that  the  presence  of  the  Ridge  Limestone  unit  causes  the  phreatic  surface  to  be 
higher  on  left  side  (right  side  of  the  plot)  of  the  valley  than  on  the  right  side 
by  10  to  IS  ft  as  was  also  noted  earlier. 
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3  Conclusions  and 
Recommendations 


The  3-D  FE  analyses  provided  a  rational  method  for  gaining  a  good  under¬ 
standing  of  the  seepage  flow  characteristics  in  the  abutments  and  foundation  of 
Cerrillos  Dam.  The  complicated  geologic  conditions  in  the  vicinity  of  the 
Dam  require  a  3-D  modeling  capability.  The  supercomputer  can  quickly 
process  data  sets  of  up  to  10,000  nodes  and  elements  (100,000  nodes  and 
elements  could  be  analyzed)  and  the  fiber  optic  network  permits  a  rapid  means 
for  transmitting  data  back  to  a  personal  computer  for  other  specialized  treat¬ 
ment.  The  adaptation  of  the  3-D  FE  program  for  other  dams  or  situations  will 
require  the  writing  of  a  grid  generator  program  (similar  to  the  MESH  pro¬ 
gram).  It  is  much  too  laborious  and  error-prone  to  attempt  to  supply  die 
location  of  the  nodes  and  elements  by  a  manual  process. 

The  use  of  3-D  FE  analyses  to  evaluate  the  effectiveness  of  proposed  cor¬ 
rective  actions  (such  as  shotcreting  or  grouting  the  permeable  units)  is  an 
effective  engineering  approach.  Two-dimensional  analysis  can  not  address 
many  issues  of  this  type. 

Because  of  the  availability  of  3-D  FE  computer  programs,  it  is  recommen¬ 
ded  that  analyses  of  these  types  be  conducted  for  all  projects  which  present  the 
possibility  of  a  3-D  effect  on  design  computations.  Even  for  cases  involving 
seepage  that  may  adequately  be  formulated  in  2-D  form,  the  use  of  numerical 
models  (2-D  FE  models)  should  receive  serious  consideration  for  use.  These 
numerical  models  generally  do  not  place  restrictions  on  the  manner  in  which 
different  permeability  zones  are  conceptually  assimilated.  Conceptual  restric¬ 
tions  generally  do  apply  for  manually  constructing  flow  nets,  especially  if  the 
transformed  section  must  be  considered.  That  is  to  say  that  flow  net  construc¬ 
tion  techniques  may  require  unnecessary  oversimplification  of  the  conceptual 
permeability  model.  For  situations  where  multiple  zones  and  anisotropic 
permeability  must  be  included,  a  numerical  model  is,  for  all  practical  pur¬ 
poses,  mandated. 
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igure  1 .  View  of  Cerrillos  Dam 
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Figure  4. 


Dam  cross  section  and  permeability  coefficients  for  2-D  flow  net  analysis 


Figure  5.  Digitized  terrain  near  Cerrillos  Dam.  (View  to 
north) 


Figure  6.  Topography  with  dam  and  location  of  Ridge 
Limestone 
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Figure  7.  Simplified  3-D  model  dam,  grout  curtain  and 
Ridge  Limestone 
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Figure  8.  View  showing  geology.  (Rockfill  shells  shown  as 
being  transparent) 


Seepage  Collection  Pipe 


Figure  9.  Close-up  view  showing  seepage  collection 
pipe 
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Figure  10.  Cut-away  view  of  dam 


Figure  1 1 .  Overhead  view  of  dam 
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Figure  12.  Section  through  dam  and  grout  curtain  midway 
up  left  valley  wall 


Figure  13.  Illustration  to  indicated  seepage  path  below  grout  curtain 


Figure  14.  3-D  finite  element  mesh 


Figure  1 5.  Pertinent  dimensions  and  features  of  3-D  FE  mesh 


klimestone  =  5  x  10-3  cm/sec 
kother  rock  =  1  x  10-4  cm/sec 
kgrout  =  1  x  10-6  cm/sec 
kdam  core  =  6.4  x  10-7  cm/sec 
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Permeable  rock  units 
extend  down  to  EL  -105  feet. 
(400  feet  below  base  of  dam) 
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Figure  1 6.  Predicted  seepage  discharges  from  3-0  FE  analysis 
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Figure  1 7.  Effect  of  placement  of  shotcrete  on  left  valley  wall 


Figure  1 8.  Effect  of  deepening  the  grout  curtain 
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Figure  19.  Piezometric  heads  from  3-D  FE 
analysis  -  pool  el  =  500  ft 


Figure  20.  Piezometric  heads  across  section 
through  dam  and  grout  curtain 


Figure  21 .  Flow  velocities  from  3-D  FE  analysis 


Figure  22.  Vertical  cross  section  showing  flow  veloc¬ 
ities  in  Ridge  Limestone  unit 


Figure  23.  Horizontal  section  near  base  of  dam 
showing  flow  velocities 


Figure  25.  Contours  of  piezometric  head  at  el  245  ft 


Figure  29.  Direction  of  seepage  flow  on  horizontal  section  at  el  370  ft 


Figure  31 .  Direction  of  seepage  flow  on  vertical  section  at  x  =  2,034  ft 


Figure  32.  Direction  of  seepage  flow  on  vertical  section  at  x  =  2,320  ft 


Figure  33.  Direction  of  seepage  flow  on  vertical  section  250  ft  downstream 
of  dam  core 


Figure  34.  Direction  of  seepage  flow  in  vertical  section  250  ft  upstream  of 
dam  core 


Figure  35.  Direction  of  seepage  flow  in  vertical  section  through  axis  of  dam 
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Figure  36.  Contours  of  phreatic  surface 


Figure  37.  Selected  cross  valley  sections  through  phreatic  surface 
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Auxiliary  Computer  Programs 
For  Cerrillos  Dam  Seepage 
Study 


MESH  -  FE  mesh  generator  program .  A2 

HEADS  -  Program  to  extract  piezometric  heads  from  FE  output  file  .  .  .  A 17 

PREFLOW  -  Program  to  extract  flow  info  from  FE  output  file . A20 

FLOWS  -  Program  to  prepare  flow  velocity  vector  files . A22 

FLOVEC  -  Program  to  draw  flow  velocity  vectors  . A25 

PHREAT  -  Program  to  extract  phreatic  surface  info  FE  output  files  ....  A32 


Appendix  A  Auxiliary  Computer  Program*  for  Carrillo*  Dam  Seepage  Study 


MESH  -  FE  mesh  generator  program 
January  1993 

PARAMETER  (NPSIZE=11560 , MSIZE=9150 ,  JSIZE=500) 

INTEGER  BLKNOD 

DOUBLE  PRECISION  XQ, YQ, ZQ, A1 , A2 , XC,  YC,  ZC 

COMMON/COOR/  X(NPSIZE) ,Y(NPSIZE)  , Z(NPSIZE)  ,  BLKNOD (MSIZE, 9) 

COMMON/ VOLINX /XQ ( 4 ) , YQ ( 4 ) , ZQ ( 4 ) 

DIMENSION  IU (NPSIZE) , IBC(JSIZE)  ,  JBC(JSIZE) ,KBC(JSIZE)  , PERM(12)  .MAT 
+TOT (12) 

CHARACTER  MATDES  (12)  *20,CODE{5)  *1 ,  CDEP*3  ,  EXT*  4  ,  HMAT*1 
CHARACTER  JASK*  1 ,  TITLE*  7 0 ,  LCDES *  7 ,  OUTFIL*  1 5 ,  GEOL*  1 ,  C*  1 ,  S ASK*  1 
CHARACTER  SHASK*1 ,  SHUP*1 ,  SHDOWN*l ,  ADUM*1 
DATA  MATDES/ 'UPSTREAM  SHELL  ' , ' SHOTCRETE  CAP 

+  ' IMPERVIOUS  CORE  ' ,  '  GROUT  CURTAIN 

+  'RIDGE  LIMESTONE  UNIT', 'UNITS  ABOVE  800  FT  ', 

+  'UNITS  BTW  750-800  FT', 'UNITS  BTW  700-800  FT', 

+  'UNITS  BTW  650-700  FT', 'UNITS  BTW  600-750  FT', 

+  'UNITS  BELOW  600  FT  ',' IMPERMEABLE  UNITS  '/ 

EXT= ' .sep' 

GO  TO  17 

14  PRINT*, 'CODE  NOT  CORRECT  _ TRY  AGAIN' 

17  CONTINUE 

PRINT* ,' ENTER  PROBLEM  CODE  (e.g.,  3LMHD100) ' 

READ  ' (5A1.A3) ' , (CODE (I) ,1=1,5) , CDEP 

IF  (CDEP.NE. ' 050'  . AND . CDEP . NE . ' 100'  . AND . CDEP . NE . ' 150'  .AND. CDEP.  NE 
+  . '200'  . 

+  AND. CDEP.NE. '250'  .AND. CDEP.NE.  '  300'  .AND. CDEP.NE. '350' )  GO  TO  14 
C=CODE (1) 

IF (C . NE . ' 2 ' . AND . C .NE . ' 3 ' )  GO  TO  14 

IF(CODE(2) .NE. 'L' ,AND.CODE(2) .NE. 'U' .AND.CODE(2) .NE. 'R' )  GOTO  14 
PRINT*, 'WANT  LAYERED  (HORIZONTALLY)  SYSTEM  ?  (Y/N) ' 

READ  ' (Al) ' , HMAT 
IF (HMAT . EQ. ' y ' )  HMAT= ' Y' 

SASK= ' N' 

PRINT*, 'WANT  TO  SHOTCRETE  THE  RIDGE  UNIT  ?  (Y/N)' 

READ  ' (Al) ' , SHASK 
IF (SHASK. EQ. ' y' )  SHASK='Y' 

SHUP= ' N' 

SHDOWN= ' N' 

IF (SHASK. EQ. 'Y' )  THEN 

PRINT*, 'ONLY  ABOVE  THE  UPSTREAM  SHELL  ?  (Y/N)' 

READ  ' (Al) ' ,  SHUP 
IF (SHUP . EQ . ' y' )  SHUP= ' Y' 

IF(SHUP.NE. ' Y' )  THEN 

PRINT*, 'ONLY  BELOW  THE  UPSTREAM  SHELL  ?  (Y/N)  ' 

READ  ' (Al) ' .SHDOWN 
I F ( SHDOWN . EQ . ' y ' )  SHDOWN=’Y' 

ENDIF 

ENDIF 

IF (CODE (3) .NE. 'L' .AND. CODE (3) .NE. 'M' .AND. CODE (3) .NE. 'H' . 
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+  AND. CODE (3) .NE. ' V' .AND. CODE (3) ,NE. 'U' ) 

IF  (CODE  (4 )  .NE.  'L'  .AND.CODEU)  .  NE  .  'M'  .  AND .  CODE  ( 4  )  .  NE  .  'H'  . 
+  AND. CODE (4) .NE. ' V' .AND. CODE (4) . NE . ' U' ) 

IF (CODE (5) .NE. 'N' .AND.CODE<5) .NE. 'S' .AND.CODE(5) . NE . 'D' ) 
IF (C . EQ. '  2 '  )  CODE ( 3 ) = ' L ' 

WRITEtOUTFIL, ' ( 5A1 , A3 , A4 ) ' )  (CODE(I) , 1=1,5) ,CDEP,EXT 
SK=1 . 22E-04 

IF(CODE{2) .EQ. 'R' )  SK=l.E-04 

IF (HMAT . EQ . 'Y' .AND. CODE (2) . NE . 'R' )  SK=1.4E-03 
PERM (6 ) =SK 

IF(CODE(3) .EQ. 'L' )  SRFL=SK 
IF (CODE (3) .EQ. 'M' )  SRFL=3.*SK 
IF(CODE(3) .EQ. 'H' )  SRFL=10.*SK 
IF (CODE (3) .EQ. 'V' )  SRFL=100.*SK 
IF(CODE(3) .EQ. 'U' )  SRFL=1000 . *SK 
PERM (5) =SRFL 

IF(CODE(4) .EQ. 'L' )  SGC=.5*SK 

IF (CODE (4) . EQ . ' M' )  SGC=.2*SK 

IF(CODE(4) .EQ. 'H' )  SGC=0.1*SK 

IF(CODE(4) .EQ. 'V' )  SGC=0.01*SK 

I F ( CODE ( 4 ) .EQ. 'U' )  SGC=0 . 000001*SK 

IF (CODE (4) .EQ. 'U' .AND. CODE (2) . EQ . 'R' )  SGC=0.001*SK 

PERM (4 ) =SGC 

PERM(l) =1000. *SK 

PERM (2) =6 .4E-07 

PERM (3) =6 .4E-07 

IF (HMAT . EQ . ' Y' )  PERM ( 12 ) =1 . E-12 
I F ( HMAT . NE . ' Y ' )  PERM (7) =1 . E-12 
I F ( C . EQ . ' 2 ' )  THEN 
ENDIF 

IF (HMAT. EQ. 'Y' )  THEN 
PERM (7)=7.8E-04 
PERM ( 8 ) =4 . 5E- 04 
PERM (9)=2.6E-04 
PERM (10) =1 . 4E- 04 
PERM (11) =1 . 4E-04 
IF (CODE (2) .EQ. 'R' )  THEN 
PERM (7) =PERM (6 ) 

PERM ( 8 ) =PERM ( 6 ) 

PERM ( 9 ) =PERM (6 ) 

PERM (10) =PERM ( 6 ) 

PERM (11) =PERM (6 ) 

ENDIF 

ENDIF 

823  PRINT*, 'THE  FOLLOWING  PERMEABILITIES  ARE  SET' 

PRINT* , ' 

KD=6 

I F ( HMAT . EQ . ' Y ' )  KD=12 
DO  828  1=1, KD 

PRINT  7023, I , MATDES ( I ) , PERM ( I ) 

828  CONTINUE 


GO 

TO 

14 

GO 

TO 

14 

GO 

TO 

14 
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7023  FORMAT (6H  K ( , 12 , 2H) - , A20 , 2H=  , 1PE8.2.8H  CM/SEC) 

PRINT*, 'TO  CHANGE  A  PERMEABILITY  VALUE  ENTER  ITS  NUMBER  (1-12)  AND 
+  ' 

PRINT* K- VALUE  (CM/SEC).  ENTER  0,0  WHEN  THROUGH' 

READ* , L, PK 

IF (L. EQ. 0 .OR . L. GT . 12 )  GO  TO  824 
PERM (L) =PK 
GO  TO  823 
824  CONTINUE 

DO  S91  1=1,12 
591  PERM (I) =PERM (I) /30 . 48 
NBORS-1 
NGRID=2 

IF (NBORS . EQ . 2 )  NGRID=1 
DATUM=850 . 

READ(CDEP, ' (13) ' )  ME LEV 
ELEV=MELEV 

PRINT* ,' HEIGHT  (FT)  OF  WATER  IN  POOL  =  ' .MELEV, '  FT' 

ELEV=ELEV+DATUM 

NELEV=ELEV 

NUMLAY= ( (NELEV+ (25*NGRID+1) ) / (50*NGRID) ) *2+1 
IF (CODE (5) .EQ. 'N' )  LC=1 
IF (CODE (5) -EQ. 'S' )  LC=2 
IF (CODE (5) .EQ. 'D' )  LC=3 

IF(LC.EQ.l)  PRINT* ,' GROUT  CURTAIN  DOES  NOT  EXIST' 

IF (LC . EQ. 2 )  PRINT* ,' GROUT  CURTAIN  IS  SHALLOW  (100  FT  DEEP)' 

IF (LC . EQ. 3)  PRINT* ,' GROUT  CURTAIN  IS  DEEP  (200  FT  DEEP)' 

PRINT* ,' OUTPUT  FILE  NAME  IS  -  ',OUTFIL 
IF (C . EQ. ' 2 ' )  PRINT*. 'THIS  IS  A  2-D  PROBLEM' 

DO  1  1=1 , NPSIZE 

1  IU(I)=0 

DO  2  1=1 , MSIZE 

2  BLKNOD ( 1 , 9 ) =  0 
NX=46/NGRID 
NY=42/NGRID 
NZ-24/NGRID 
XDEL=  71.4  *NGRID 
YDEL=50.*NGRID 
ZDEL=25 . *NGRID 
NP=1 

DO  100  1=1, NUMLAY 

MY=NY 

YADD=  0 . 

IF(I- ( (I/2)*2) . EQ . 0)  THEN 

YADD=2S.*NGRID 

MY-NY-1 

END  IF 

DO  200  J=1 , MY 

MX=NX 

XADD-0. 

IF(I- ( (1/2) *2) .EQ.O)  THEN 
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XADD=35.7*NGRID 

MX=NX-1 

ENDIF 

DO  300  K=1 , MX 
Z ( NP )  =  ( I  - 1 )  *  ZDEL 
Y(NP) = (J-l) *YDEL+YADD 
X (NP) = (K-l) *XDEL+XADD 
NP=NP+1 
300  CONTINUE 
200  CONTINUE 
100  CONTINUE 
NP=NP-1 

DO  588  I=NP+1 , NPSIZE 
588  IU(I)=1 
50  FORMAT (2IS,4F10. 2) 

NELX = NX- 3 
NELY=NY-2 
NEL=1 

DO  400  K=1 , NUMLAY/2 
NL=  (K-l)  *  (  (NX*NY)  +  (  (NX-1)  *  (NY-1)  )  ) 
NU=NL+NX*NY 
DO  500  1=1 , NELY 
DO  600  J=1 , NELX 
BLKNOD ( NEL , 1 ) =NL+J+1+ (1-1)  *NX 
BLKNOD(NEL, 2) =NL+J+2+ (1-1)  *NX 
BLKNOD (NEL, 3) =NL+J+3+I*NX 
BLKNOD (NEL, 4) =NL+J+2+I*NX 
BLKNOD (NEL, 5) =NU+J+ (1-1) * (NX-1) 
BLKNOD (NEL, 6) =NU+J+1+ ( I - 1 ) * (NX-1) 
BLKNOD (NEL, 7) =NU+J+2+I* (NX-1) 
BLKNOD (NEL, 8) =NU+J+1+I* (NX-1) 
NEL=NEL+1 
600  CONTINUE 
500  CONTINUE 
NL=NU 

NU=NL+ (NX-1) * (NY-1) 

DO  550  1=1, NELY 
DO  650  J=1 , NELX 

BLKNOD (NEL, 1) =NL+J+ (1-1) * (NX-1) 
BLKNOD (NEL , 2 ) =NL+J+1+ ( I- 1) * (NX-1) 
BLKNOD (NEL, 3) -NL+J+2+I* (NX-1) 
BLKNOD (NEL, 4) -NL+J+l+I* (NX-1) 
BLKNOD (NEL, 5) =NU+J+I*NX 
BLKNOD (NEL, 6) =NU+J+1+I*NX 
BLKNOD (NEL, 7) =NU+J+2+ (1+1) *NX 
BLKNOD (NEL, 8) =NU+J+1+ (I+I) *NX 
NEL=NEL+1 
650  CONTINUE 
550  CONTINUE 
400  CONTINUE 
NEL=NEL- 1 
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7  FORMAT (15, 3F15 . 1) 

PRINT* , ' #  NODES = ' , NP , '  #  ELEMENTS = ' , NEL 

C- COMPUTE  BANDWIDTH  OF  ORIGINAL  ELEMENT  MESH 
MAX=0 

DO  7345  1=1, NEL 

MAX=MAXO (BLKNOD (1,7) -BLKNOD (1,1) , MAX) 

7345  CONTINUE 
C-DO  GEOLOGY  LAYERS 
GEOL= ' Y ' 

IF(GEOL.EQ. 'y' . OR . GEOL . EQ . ' Y ' )  THEN 

PRINT*, 'WAIT . ASSIGNING  MATERIAL  NUMBERS  TO  GEOLOGIC  LAYERS' 

IFIND=0 
LCOUNT=0 
NCOUNT=0 
DX=-71.4*NGRID 
DY-50 . *NGRID 
DZ=- 50 . *NGRID 
XD=3141 . 6 
YD=0  . 

ZD=1200 . 

DO  77  J=l, 500 

IF (NBORS . EQ . 1 . AND . ( J+3- (NGRID-1) ) . LT . 22 )  GO  TO  77 

IF (NBORS . EQ . 1 . AND  .  (J+3 - (NGRID- 1) )  . GT . 24 )  GO  TO  79 

IFtNBORS. EQ. 2. AND. (J+3- (NGRID-1) ) .LT. 44)  GO  TO  77 

IF(NBORS. EQ. 2. AND. (J+3- (NGRID-1) ) .GT. 49)  GOTO  79 

NCOUNT=0 

DO  78  1=1, NEL 

II*I 

CALL  CENTER (II, XQ ( 4 ) , YQ ( 4 ) , ZQ ( 4 ) ) 

XQ (2 ) =XD+ (J-l) *DX 
XQ(1) =3141.6 
XQ ( 3 ) =XQ ( 1 ) 

YQ ( 1 ) =YD+ (J-l) *DY 
YQ(2) =0. 

YQ(3) =YQ(2) 

ZQ(3)=ZD+ (J-l) *DZ 
ZQ(1) =1200. 

ZQ  ( 2 )  =  ZQ  ( 1 ) 

CALL  VOL(Al) 

IF(J.EQ.l)  Al=l . 

IF(I.EQ.O)  THEN 
PRINT*, ' VOL= ' , A1 

PRINT  9 9 , XQ ( 1 ) , XQ ( 2 ) ,XQ(3) ,YQ(1) ,YQ(2) ,YQ(3)  ,ZQ(1)  ,ZQ(2)  ,ZQ(3) 
ENDIF 

99  FORMAT (4H  XQ=,9F7.1) 

XQ (2) =XQ (2) +DX 
YQ ( 1) =YQ (1) +DY 
ZQ(3) =ZQ(3) +DZ 
CALL  VOL ( A2 ) 

IF(I.EQ.O)  THEN 
PRINT*, ' VOL=' , A2 
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PRINT  9 9 , XQ ( 1 ) , XQ ( 2 )  , XQ (3)  ,YQ(1)  ,YQ(2)  ,YQ(3) , ZQ ( 1 ) ,ZQ(2)  , ZQ ( 3 ) 
ENDIF 

IF (A1*A2 . LT . 0 . )  THEN 
IF (IFIND.EQ. 0)  JJ=J 
IFIND=1 

BLKNOD (1,9) =J+3 - (NGRID-1) 

NCOUNT=NCOUNT+l 
LCOUNT=LCOUNT+ 1 
ENDIF 

78  CONTINUE 
77  CONTINUE 

79  CONTINUE 
ENDIF 

PRINT* ,  ' REDUCE  GEOLOGIC  LAYERS  TO  TWO  MATERIALS' 

PRINT* MATERIAL  #  5  IS  THE  RIDGE  LIMESTONE' 

PRINT* MATERIAL  #  6  IS  ASSIGNED  TO  ALL  OTHER  STRATA' 

DO  8834  1  =  1 , NEL 
MMM=BLKNOD (1,9) 

BLKNOD (I, 9) =6 

IF (NBORS . EQ . 1 . AND . MMM . GE . 22 . AND . MMM. LE . 24 )  BLKNOD (I, 9) =5 
IF (NBORS . EQ . 2 . AND . MMM . GE . 44 .AND. MMM. LE. 49)  BLKNOD (I, 9) =5 
8834  CONTINUE 

C-SET  MATERIAL  NUMBER  TO  1000  IN  ALL  VALLEY  ELEMENTS 
YQ ( 1 ) =0. 

ZQ(1) =850. 

YQ (2 ) =0 . 

ZQ(2) =1200. 

YQ(3) =100. 

ZQ(3) =1200. 

NCOUNT=0 
LCOUNT  =  0 
DO  2200  1=1, NEL 
XQ (1) =1785. 

XQ (2) =2284.8 
XQ ( 3 ) =2284 . 8 

CALL  CENTER ( I , XQ ( 4 ) , YQ ( 4 ) , ZQ ( 4 ) ) 

CALL  VOL(Al) 

XQU)  =135^5 .6 
XQ(2) =856.8 
XQ(3) =856 .8 
CALL  VOL ( A2 ) 

IF  (A1*A2 . LT .0.0. AND . ZQ (4 )  .GT.850.)  THEN 
BLKNOD (I, 9) =1000 

IF (YQ (4 ) .GT. 900.0. AND. YQ(4) .LT.1050.)  THEN 

BLKNOD (I, 9) =3 

LCOUNT= LCOUNT+ 1 

GO  TO  2200 

ENDIF 

NCOUNT =NCOUNT + 1 
ENDIF 

2200  CONTINUE 
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MCOUNT=NCOUNT 

PRINT*, 'ASSIGNED' ,  LCOUNT,  '  ELEMENTS  IN  DAM  CORE  TO . MATER  I 

+AL  #3' 

IF (LC.GT. 1)  THEN 
LCOUNT =0 
NCOUNT=0 
YQ(1)=0. 

ZQ  (1) *650 . 

IF (LC . EQ. 2 )  ZQ ( 1 ) =750. 

YQ(2)=100. 

ZQ{2) =650. 

IF (LC. EQ. 2 )  ZQ (2 ) =750 . 

YQ(3)=0. 

ZQ (3) =1200 . 

DO  66  1=1, NEL 

CALL  CENTER (I, XQ (4) ,YQ(4) ,ZQ(4)  ) 

XQ (1) =1785 . 

XQ(2) =1785. 

XQ (3) =2570 . 4 

IF (LC.EQ. 2)  XQ (3 ) =2427 . 6 

CALL  VOL(Al) 

XQ (1) =1356 . 6 
XQ (2) =1356 . 6 
XQ(3) =571.2 

IF (LC.EQ. 2)  XQ (3) =714.0 
CALL  VOL ( A2 ) 

IF  (A1*A2 . LT .0.0 . AND . ZQ (4 }  .GT.ZQ(l) )  THEN 

IF  (ABS (BLKNOD (1,9) )  .NE . 3 )  THEN 

IF (YQ (4 ) .GT. 900.0. AND. YQ(4) .LT.1050.)  THEN 

BLKNOD (I, 9) =4 

NCOUNT=NCOUNT+ 1 

ENDIF 

ENDIF 

ENDIF 

66  CONTINUE 

PRINT*, 'ASSIGNED' ,  NCOUNT,  '  ELEMENTS  IN  GROUT  CURTAIN  TO....MATERI 
+AL  #4' 

ENDIF 
NCOUNT=0 
LCOUNT =0 
XQ (1) =856 . 8 
XQ (2 ) =2284 . 8 
XQ(3) =1356 .6 
ZQ (1) =1200 . 

ZQ (2 ) =1200 . 

ZQ(3) =850. 

DO  366  1=1, NEL 

IF (BLKNOD (I, 9) .EQ. 1000)  THEN 

YQ (1) =800. 

YQ (2) =800 . 

YQ(3) =500. 
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CALL  CENTER (I, XQ (4) ,YQ(4) , ZQ(4) ) 

CALL  VOL(Al) 

YQ(1) =1100. 

YQ (2 ) =1100. 

YQ(3) =1400. 

CALL  VOL ( A2 ) 

IF (A1*A2 .LT. 0 . )  THEN 

ENDIF 

ENDIF 

366  CONTINUE 

PRINT* FOUND  ' , MCOUNT - NCOUNT - LCOUNT , '  ELEMENTS  IN  VALLEY.' 

IF (C . EQ . ' 2 ' )  THEN 

NCOUNT=Q 

DO  9234  1  =  1 , NEL 

CALL  CENTER (I, XC, YC, ZC) 

IF (XC . LT. 1356.6. OR . XC . GT . 1785 . )  THEN 
IF (HMAT.NE. ' Y' )  BLKNOD ( I , 9) =7 
IF(HMAT.EQ. ' Y' )  BLKNOD (I , 9) =12 
NCOUNT=NCOUNT+ 1 
ENDIF 

9234  CONTINUE 

PRINT* , NCOUNT, '  ELEMENTS  OUTSIDE  2-D  ZONE  GIVEN  VERY  HIGH  PERMEABI 
+LITY . ' 

ENDIF 

I F ( HMAT . EQ . ' Y ' )  THEN 
DO  9898  1=1, NEL 
C=CODE (1) 

CALL  CENTER ( I, XC.YC.ZC) 

IF (ZC.LT. 600. AND. BLKNOD (I, 9) . LE . 6 . AND . C . EQ. ' 2 ' )  BLKNOD (I, 9) =11 
IF (ZC.LT. 600. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 3' )  BLKNOD (I, 9) =11 
IF (ZC.LT. 6 50. AND. ZC.GT. 600. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 3' )  BLKN 
+OD (1,9) =10 

IF (ZC.LT. 6 50 .AND. ZC.GT. 600. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 2' )  BLKN 
+OD (1,9) =10 

IF (ZC.LT. 6 50. AND. ZC.GT. 600. AND. BLKNOD (I, 9) . EQ. 5 . Alto. C. EQ . '2' )  BLKN 
+OD (1,9) =10 

IF (ZC.LT. 700. AND. ZC.GT. 6 50. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . '3' )  BLKN 
+OD (1,9) =9 

IF (ZC.LT. 700. AND. ZC.GT. 6 50. AND. BLKNOD (I, 9) .EQ. 6 .AND.C.EQ. '2' )  BLKN 
+OD (1,9) =9 

IF (ZC.LT. 700. AND. ZC.GT. 6 50. AND. BLKNOD (I, 9) . EQ . 5 . AND . C . EQ . '2')  BLKN 
+OD (1,9) =9 

IF (ZC.LT. 750. AND. ZC.GT. 7 00. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 3' )  BLKN 
+OD (I , 9) =8 

IF (ZC.LT. 750. AND. ZC.GT. 7 00. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 2' )  BLKN 
+OD (1,9) =8 

IF (ZC.LT. 750. AND. ZC.GT. 7 00. AND. BLKNOD (I, 9) . EQ . 5 . AND . C . EQ . '2' )  BLKN 
+OD ( I , 9) =8 

IF (ZC.LT. 800. AND. ZC.GT. 7 50. AND. BLKNOD (I, 9) . EQ . 6 . AND . C . EQ . ' 3' )  BLKN 
+OD (1,9) =7 
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IF (ZC.LT. 800 . AND. ZC.GT. 750. AND. BLKNOD(I, 9)  . EQ.  6 . AND.C. EQ.  '  2  '  )  BLKN 
+OD { I , 9) =7 

IF  (ZC.  LT.  800  .  AND .  ZC  .GT  .750  .  AND .  BLKNOD  ( 1 , 9 )  .  EQ .  5  .  AND .  C  .  EQ .  '2'  )  BLKN 
+OD(I,9)=7 

IF (ZC . LT. 450 . AND . CODE (2 ) . EQ . ' R ' )  BLKNOD ( I , 9) =12 
9898  CONTINUE 
END  IF 

IF (SASK.EQ . ' L' )  THEN 

DO  6490  I«1 , NEL 

CALL  CENTER { I , XC , YC , ZC ) 

IF (XC.GT. 1S70 . 8)  BLKNOD (1 , 9)  =12 

6490  CONTINUE 
END  IF 

IF (SASK.EQ. ' R' )  THEN 

DO  6491  1=1, NEL 

CALL  CENTER ( I , XC , YC , ZC ) 

IF (XC.LT. 1S70 . 8)  BLKNOD (I , 9) =12 

6491  CONTINUE 
END  IF 

IF (2 . EQ. 2}  THEN 

C-WRITE  ORIGINAL  NODE  USE  FILE 
DO  478  1=1, NEL 
DO  478  K=1 , 8 
KKK=BLKNOD ( I , K) 

IU(KKK) =IU (KKK) +1 

478  CONTINUE 

OPEN (19, FILE= '  SCR '  ) 

DO  479  1=1, NP 

479  WRITE (19, *)  IU(I) 

C- STRIP  ALL  VALLEY  ELEMENTS  AND  ALL  ELEMENTS  ABOVE  POOL 
PRINT*, 'WAIT. . .NOW  STRIPPING  VALLEY  ELEMENTS  ' 

M=1 

990  CONTINUE 

DO  770  I=M, NEL 

CALL  CENTER  ( I , XQ (4 )  , YQ (4 )  , ZQ (4 ) ) 

IF (ABS (BLKNOD (1,9))  . EQ . 1000 .  )  THEN 
M=I 

DO  775  J=I , NEL- 1 
DO  775  K=1 , 9 

775  BLKNOD (J,K) -BLKNOD (J+l.K) 

NEL-NEL- 1 
GO  TO  990 
ENDIF 

770  CONTINUE 

PRINT* ,  '  TOTAL  NUMBER  OF  ELEMENTS  NOW  '  ,  NEL 
2 -COUNT  NUMBER  OF  ELEMENTS  IN  EACH  TYPE  OF  MATERIAL 
DO  8601  1=1,12 
8601  MATTOT (I) =0 

DO  8602  1=1, NEL 
KK= BLKNOD (I, 9) 

MATTOT ( KK) -MATTOT ( KK) + 1 
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8602  CONTINUE 

DO  8603  1  =  1, 12 

8603  PRINT  8604 , MATTOT (I) , I 

8604  FORMAT ( 110, 32H  ELEMENTS  ASSIGNED  TO  MATERIAL  #,I3) 

DO  375  1=1 , NPSIZE 

375  IU(I)=0 

DO  88  1=1, NEL 
DO  88  K=1 , 8 
KKK=BLKNOD ( I , K) 

IU(KKK) =1 

88  CONTINUE 

OPEN ( 18 , FILE= ' NOTUSE' ) 

FEWIND  18 

WRITE (18,*)  'LIST  OF  UNUSED  NODES' 

DO  89  1=1, NP 

IF ( IU ( I) . EQ . 0)  WRITE (18 , 7)  I,X(I) , Y (I) , Z (I) 

89  CONTINUE 
CLOSE  (18) 

C-STRIP  UNUSED  NODES  FROM  FILE  'SCR' 

REWIND  19 

OPEN (22, FILE=' IUS ' ) 

M=0 

DO  913  1=1 , NP 
READ (19,*)  JJ 
IF(IU(I) .GT.0)  THEN 
M=M+1 

WRITE (22,*)  JJ 
ENDIF 

913  CONTINUE 
CLOSE (22) 

CLOSE (19) 

C-STRIP  UNUSED  NODES  AND  RENUMBER  NODES  AND  ELEMENTS 

PRINT* , ' WAIT . . . REMOVING  UNUSED  NODES  AND  RENUMBERING  ELEMENTS ' 
NPP=NP 
IR  =  0 
11  =  0 

DO  2880  1=1, NPP 
IF(IU(I) .EQ. 0)  THEN 
11=11+1 

DO  2882  K=I - IR , NP- 1 
X (K) =X (K+l ) 

Y (K) =Y (K+l ) 

Z (K) =Z (K+l) 

2882  CONTINUE 

DO  2883  L=-l ,  NEL 

IF (BLKNOD (L, 7) .LT.I-IR)  GO  TO  2883 
DO  2884  K= 1 , 8 

IF (BLKNOD (L,K)  .GE.I-IR)  BLKNOD (L, K) =BLKNOD (L, K) -1 
2884  CONTINUE 

2883  CONTINUE 
IR= IR+ 1 
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NP=NP- 1 
ENDIF 

2880  CONTINUE 

PRINT* , ' TOTAL  NUMBER  OF  NODES  NOW  ' , NP 
ENDIF 

DO  444  1=1, NP 
444  IUCI)=0 

DO  188  1=1, NEL 
DO  188  K=1 , 8 
KKK=BLKNOD ( I , K) 

IF {KKK . LT . 1 . OR . KKK. GT .NPSIZE)  PRINT* , KKK 
IU  (KICK)  =1 

188  CONTINUE 

OPEN ( 18 , FILE= ' NOUSE ' ) 

REWIND  18 

WRITE (18,*)  'LIST  OF  UNUSED  NODES' 

DO  189  1=1, NP 

IF(IU(I).EQ.O)  WRITE (18,7)  I, X (I) ,Y(I) ,Z(I) 

189  CONTINUE 
CLOSE (18) 

C- COUNT  NUMBER  OF  TIMES  EACH  NODE  USED  TO  FORM  AN  ELEMENT 
DO  544  1=1, NP 
544  IU(I)=0 

DO  538  1=1, NEL 
DO  538  K=1 , 8 
KKK=BLKNOD ( I , K) 

IU(KKK) =IU (KKK) +1 
538  CONTINUE 

PRINT* ,  ' WAIT ....  LOCATING  NODES  TO  APPLY  NODAL  PRESSURES  ' 
C-FIND  NODES  ON  UPSTREAM  SIDE  OF  DAM  TO  APPLY  NODAL  PRESSURES 
OPEN (22 , FILE= ' IUS ' ) 

REWIND  22 
NUMBI=0 
NUMBJ=0 
NUMBK=0 

DO  3000  1=1, NP 
READ (22,*)  IUJ 
IF(IUJ.NE. IU(I) )  THEN 

IF(Yd)  .GT.  950 . 0  .  AND.  Z  (I)  .LE.ELEV+1.)  THEN 
NUMB I =NUMB I + 1 
IBC(NUMBI) =1 
ENDIF 

IF (Y(I) .LT. 950.0. AND. Z(I) .LE.851.)  THEN 
NUMB  J = NUMB J + 1 
JBC (NUMBJ) =1 
ENDIF 

IF (Y (I ) .LT. 950.0. AND. Z(I) .GT.851.)  THEN 

NUMBK=NUMBK+1 

KBC(NUMBK) =1 

ENDIF 

ENDIF 
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3000  CONTINUE 
CLOSE (22) 

MAX=0 

DO  7692  I=1,NEL 

MAX = MAX 0 (BLKNOD ( 1 , 7) -BLKNOD(I.l)  .MAX) 

7692  CONTINUE 

PRINT*,  'WAIT.  ..  .  WRITING  NODE  AND  ELEMENT  DATA  TO  FILE  ’.OUTFIL 
OPEN (15, FILE='UPNODE'  ) 

DO  778  1=1 , NUMBI 

WRITE (15, 63)  IBC(I)  ,X(IBC(I) )  ,Y(IBC(I) ) , Z(IBC(I) )  ,  ELEV 
63  FORMAT (16, 3 FI 0. 2, F2 0.1) 

778  CONTINUE 
FLEV=850 . 

DO  779  1=1 , NUMBJ 

WRITE (15,63)  JBC(I) ,X(JBC(I) ) ,Y(JBC(I) ) ,Z(JBC(I) ) , FLEV 

779  CONTINUE 

DO  781  1=1 , NUMBK 

WRITE (15, 63)  KBC(I)  ,X(KBC(I)  )  ,  Y(KBC(I) ) , Z(KBC(I)  )  , Z (KBC(I) ) 

781  CONTINUE 
CLOSE (15) 

OPEN ( 14 , FILE= ' USED ' ) 

WRITE (14,58)  (I, IU(I) , 1=1, NP) 

58  FORMAT (5 (18, 1H ( , II, 1H) ) ) 

CLOSE (14) 

GO  TO  10001 
10002  CONTINUE 

DO  1800  1=1, NP 
IU(I) =0 
BDVL=0 . 0 

DO  1801  J=l, NUMBI 
IF ( IBC ( J) . EQ. I )  THEN 
IU(I) =1 

BDVL= ELEV - DATUM 
GO  TO  1802 
ENDIF 

1801  CONTINUE 

DO  1803  K=l, NUMBJ 
IF(JBC(K) .EQ.I)  THEN 
IU ( I ) =1 
BDVL=0 . 

GO  TO  1802 
ENDIF 

1803  CONTINUE 

DO  1805  L=l, NUMBK 
IF(KBC(L) .EQ.I)  THEN 
IU(I) =2 

BDVL=2 (I) -DATUM 
GO  TO  1802 
ENDIF 

1805  CONTINUE 

1802  CONTINUE 
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WRITE (13 ,50)  I , IU (I) ,X(I) ,Y(I) ,Z(I) , BDVL 

1800  CONTINUE 

C  SET  MATERIAL  NUMBER  OF  CAPPED  ELEMENTS 
I COUNT =0 

IF(SHASK.EQ. 'V' )  THEN 
OPEN (55 , FILE= ' BOSSCAP' ) 

PRINT* , ' . - . WAIT  FINDING  SHOTCRETED  ELEMENTS ' 

DO  109  1=6000 , NEL 

IF (BLKNOD (I , 9) .EQ.5)  THEN 

IF (Z (BLKNOD (1,8)) . LT .851. )  GO  TO  109 

DO  110  K=1 , 8 

N=BLKNOD { I , K) 

DO  110  J=1 , NP 

IF(Z(J) .LT.851.)  GO  TO  110 

IF(N.EQ.J)  THEN 

IF (IU (J) . EQ. 1)  THEN 

Al-1. 

A2  =  - 1 . 

IF (SHUP . EQ. 'Y'  . OR . SHDOWN . EQ .  '  Y'  )  THEN 
CALL  CENTER ( I , XQ ( 4 ) , YQ ( 4 ) , ZQ ( 4 ) ) 

ZQ (1) =1200. 

ZQ (2 ) =1200 . 

ZQ(3) =850. 

XQ (1 ) =856 . 8 
XQ (2) =2284 . 8 
XQ (3) =1356 . 8 
YQ ( 1 ) =1100. 

YQ (2 ) =1100 . 

YQ{3) =1400. 

CALL  VOL(Al) 

YQ(1) =3000. 

YQ(2) =3000. 

YQ(3) =3000. 

CALL  VOL ( A2 ) 

ENDIF 

IF  (SHDOWN. NE  .  '  Y'  .  AND  .  SHUP  .  NE  .  'Y'  )  THEN 
BLKNOD (I, 9) =2 
I COUNT = I COUNT +1 
XX=XQ (4 ) 

YY=YQ (4 ) 

ZZ=ZQ (4 ) 

WRITE ( 55 , * )  I , XX , YY , ZZ 

GO  TO  109 

ENDIF 

IF (SHDOWN. EQ . ' Y'  . AND. A1*A2 .GE  .  0  .  )  THEN 
BLKNOD (I, 9) =2 
I COUNT = I COUNT + 1 
XX=XQ (4 ) 

YY=YQ (4 ) 

ZZ=ZQ (4 ) 

WRITE ( 55 , * )  I , XX , YY , ZZ 


A14 


Appendix  A  Auxiliary  Computer  Program*  for  Carrillo*  Dam  Seepage  Study 


GO  TO  109 
ENDIF 

IF (SHUP . EQ . ' Y' . AND .  A1  *  A2  .  LT . 0 . )  THEN 
BLKNOD (1,9) =  2 
I  COUNT  =  I  COUNT +1 
XX=XQ(4) 

YY=YQ (4 ) 

ZZ*ZQ(4) 

WRITE (55 , *)  I ,  XX, YY, ZZ 

GO  TO  109 

ENDIF 

ENDIF 

ENDIF 

110  CONTINUE 
ENDIF 

109  CONTINUE 
CLOSE (55) 

PRINT*, ICOUNT, '  ELEMENTS  ASSIGNED  TO  MATERIAL  #2' 

ENDIF 

DO  2000  1*1, NEL 

2000  WRITE (13,55)  I ,  (BLKNOD ( I , K)  ,  K=1 , 9 ) 

55  FORMAT (1015) 

GO  TO  10003 
10001  CONTINUE 
MATNO* 0 
DO  378  1*1 , NEL 

IF (BLKNOD (1,9)  . GT . MATNO)  MATNO*  BLKNOD (1,9) 

378  CONTINUE 

OPEN (13 , FILE-OUTFIL) 

IF(LC.EQ.l)  LCDES* '  NO' 

I F ( LC . EQ . 2 )  LCDES* 'SHALLOW' 

IF (LC . EQ . 3 )  LCDES*'  DEEP' 

NOVEL* 0 

WRITE (13,315)  ELEV-DATUM , LCDES  ,  OUTFIL 
WRITE ( 13 , ' (4I5.F10.2) ' )  NP , NEL , MATNO , NOVEL , DATUM 
315  FORMAT (16H  HEIGHT  OF  POOL* , F10 . 0 , 3X,  A7 , 14H  GROUT  CURTAIN, 6H  FILE  , 
+A15 ) 

XK=1 . 0 
YK=1 . 0 
ZK=1 . 0 

DO  918  1*1, MATNO 

918  WRITE (13, ' ( 15 , 3E10 . 3 ) ' )  I , PERM ( I ) , PERM ( I ) , PERM ( I ) 

GO  TO  10002 
10003  CONTINUE 
CLOSE (13) 

C-REMOVE  TEMP  FILES 

OPEN (13, FILE* 'SCR'  ) 

CLOSE (13, STATUS* 'DELETE'  ) 

OPEN (13, FILE*' NOTUSE'  ) 

CLOSE ( 13 , STATUS* ' DELETE '  ) 

OPEN (13, FILE* ' IUS ' ) 
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CLOSE ( 13, STATUS =' DELETE'  ) 

OPEN ( 13 , FILE= ' NOUSE ' ) 

CLOSE < 13 .STATUS-' DELETE' ) 

OPEN (13 , FILE= ' UPNODE' ) 

CLOSE { 13 , STATUS- ' DELETE '  ) 

OPEN (13 , FILE- ' USED' ) 

CLOSE { 13 , STATUS- ' DELETE '  ) 

STOP 

END 

SUBROUTINE  CENTER (I , XC, YC( ZC) 

PARAMETER  (NPSIZE-11560 , MSIZE-9150 , JSIZE-500) 

INTEGER  BLKNOD 

DOUBLE  PRECISION  XC,YC,ZC 

CO MMON/COOR/  X(NPSIZE) , Y(NPSIZE) , Z (NPSIZE) , BLKNOD (MSI ZE, 9) 

XX -0 . 

YY-0 . 

ZZ-0. 

DO  100  <J-1 , 8 
XX-XX+X (BLKNOD (I , J) ) 

YY-YY+Y (BLKNOD ( I , J) ) 

100  ZZ-ZZ+Z (BLKNOD (I, J) ) 

XC-XX/8. 

YC-YY/8 . 

ZC-ZZ/8. 

RETURN 

END 

SUBROUTINE  VOL(VL) 

DOUBLE  PRECISION  X. Y. Z, A, B, C, D, VL, XI , X2 , X3 , Y1 . Y2 , Y3 , Z1 , Z2 , Z3 
COMMON/VOLINX/X (4 ) ,Y(4)  ,Z(4) 

X1=X(1) -X (4 ) 

X2=X (2) -X (4 ) 

X3-X ( 3 ) -X (4 ) 

Y1=Y (1) -Y (4 ) 

Y2=Y (2) - Y (4 ) 

Y3=Y (3 ) - Y (4 ) 

Z1=Z(1) - Z (4 ) 

Z2=Z (2) - Z (4 ) 

Z3=Z ( 3 ) -Z (4 ) 

VL-X1* (Y2*Z3-Y3*Z2) -X2* (Y1*Z3-Y3*Z1) +X3* (Y1*Z2-Y2*Z1) 

RETURN 

END 
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$ FREE FORM 

PARAMETER  (NPSIZE-11560,NEL$IZ«9150) 

DIMENSION  X(NPSIZE) ,Y<NPSIZE) ,Z(NPSIZE> ,IBC(NPSIZE) , HEAD (NPSIZE) . - 
FLOW (NPSIZE) 

CHARACTER  FIL*12,OUTFIL(100) *5 , DUM*2 , TITLE* 76, B*2, 0*2, A*2, ZFL*1, XFL*1 ,  YFL*1 

CHARACTER  OUT(5)*l 

ZFL-'Z' 

XFL-'X' 

YFL-'Y' 

PRINT* , ' ENTER  INPUT  FILE  NAME' 

READ  ' (A12) ' , FIL 
OPEN (10 , FILE*FIL) 

CALL  REDUM(S) 

READ (10, ' (A76) ' )  TITLE 
CALL  REDUM(S) 

12  FORMAT  (SIX, 15) 

READ(10, 12)  NP 
CALL  REDUM(l) 

READ (10 , 12)  NEL 
CALL  REDUM(l) 

READ (10, 12)  MATNUM 
CALL  REDUM(l) 

READ (10, ' (SIX, Fll . 0) ' )  DATUM 
CALL  REDUM(IO) 

CALL  REDUM (MATNUM) 

CALL  REDUM (10) 

DO  100  1*1 , NP 

READ(10,*)  J,IBC(I) ,X(I) ,Y(I) ,Z(I) 

100  CONTINUE 
CALL  REDUM (NEL+10) 

CALL  REDUM (10) 

DO  200  1*1, NP 

READ(10. ' ( 16 , 2F13 . 0 , IX , 3 A2 ) ')  J,HEAD(I) ,FLOW(I) ,B,0,A 

IF (B . EQ. '  *')  IBC (I) *1 

IF (O . EQ. '  *')  IBC ( I) *2 

IF(A.EQ. '  *')  IBC (I) *3 

200  CONTINUE 

CALL  REDUM (7) 

READ (10, ' (4 OX, F12 . 0) ' )  AVFLOW 
"  WRITE  THE  Z-PLANES  FILES 
ZMAX--99999. 

XMAX-- 99999 . 

YMAX--99999. 

ZMIN--ZMAX 

XMIN--XMAX 

YMIN.-YMAX 

DO  300  I-l.NP 

I F ( Z  ( I )  . GT.ZMAX)  ZMAX»Z(I) 
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IF (X { I ) .GT.XMAX)  XMAX*X(I) 

IF(Y(I) .GT.YMAX)  YMAX=Y(I) 

IF(2(I)  .LT.ZMIN)  ZMIN=-2  ( I ) 

IF(X(I) .LT.XMIN)  XMIN=X ( I ) 

IF(Y(I) . LT.YMIN)  YMIN*Y(I) 

300  CONTINUE 
PRINT*, 'X' , XMIN, XMAX 
PRINT*, *  Y' , YMIN, YMAX 
PRINT* , ' Z ' , ZMIN , ZMAX 
J= (1 . 00l*ZMAX/50) +1 
DO  400  1=1, J 
E= (1-1) *50. 

IE= (E+ . 001) 

WRITE (OUTFIL (I) , ' (A1 , 14 ) ' )  ZFL, IE 
READ (OUTFIL ( I ) , ' (5A1) ' )  (OUT(L) ,L=1,S) 

DO  4  L=l, 5 

IF (OUT (L)  . EQ. '  '  )  OUT (L) = ' 0 ' 

4  CONTINUE 

WRITE  (OUTFIL  (I)  ,  '  (5A1)  '  )  ('  T'(L)  ,L«1,5) 

OPEN ( 12 , FILE.OUTFIL (I) ) 

DO  500  K=1 ,  NP 

IF(X(K) .GT. XMIN. AND. X(K) .LT. XMAX)  THEN 
IF(Y(K) .GT. YMIN. AND. Y(K) -LT. YMAX)  THEN 
IF (Z (K) . GT. E-l . 0 .AND. Z (K) .LT.E+1.)  THEN 

WRITE (12 , *  ( 15 , 4E12 . 4 , 15) * )  K,X(K) ,Y(K)  , HEAD (K)  ,FLOW(K)  ,IBC(K) 
END  IF 
END  IF 
ENDIF 

500  CONTINUE 
CLOSE (12) 

400  CONTINUE 
"  WRITE  Y- PLANES 
J= (1 . 001* YMAX/ 50) +1 
DO  410  1*1, J 
E= (I- 1 ) *50. 

IE= (E+0 . 001) 

WRITE (OUTFIL (I) , ’ (Al, 14) ' )  YFL, IE 
READ (OUTFIL (I) , ' (5A1) ' )  (OUT(L) ,L=1, 5) 

DO  6  L«l, 5 

IF(OUT(L) .EQ. '  ')  OUT (L) = ' 0 ' 

6  CONTINUE 

WRITE (OUTFIL (I) , ' (5A1) ' )  (OUT(L) ,L=1, 5) 

OPEN (12 , FILE=OUTFIL (I)  ) 

DO  510  K*1 ,  NP 

IF (X (K) . GT. XMIN. AND. X(K) .LT.XMAX)  THEN 
IF (Z (K) .GT. ZMIN. AND. Z (K) .LT.ZMAX)  THEN 
IF (Y (K) . GT . E- 1 . 0 . AND . Y (K) .LT.E+1.)  THEN 

WRITE (12, ' (I5.4E12.4, 15) ' )  K,X(K) ,Z(K) , HEAD (K) ,FLOW(K) ,IBC(K) 

ENDIF 

ENDIF 

ENDIF 
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510  CONTINUE 
CLOSE (12) 

410  CONTINUE 
"WRITE  X- PLANES 

(1 . 001*XMAX/142 . 8 ) +1 
DO  420  1*1,  J 
E* (1-1) *142 . 8 
IE=E+0 . 001 

WRITE (OUTFIL ( I ) , ' ( A1 , 14 ) ' )  XFL, IE 
READ (OUTFIL ( I ) , ' (SA1) ')  (OUT(L) ,L=1,5) 

DO  8  L=1 , 5 

IF (OUT (L) .  EQ.  '  ')  OUT (L) = ' 0 ' 

8  CONTINUE 

WRITE (OUTFIL (I)  ,  '  (5A1) ' )  (OUT(L) ,L=1, 5) 

OPEN ( 12 , FILE*OUTFIL (I) ) 

DO  520  K*1,NP 

IF (Y (K) . GT.YMIN.AND. Y (K) .LT.YMAX)  THEN 
IF (Z (K) . GT . ZMIN . AND . Z ( K) .LT.ZMAX)  THEN 
IF (X (K) . GT . E-l . 0 .AND. X (K) .LT.E+1.)  THEN 

WRITE (12 , ' ( 15 , 4E12 . 4 , 15) ')  K, Y (K) , Z (K) , HEAD (K) ,FLOW(K) ,IBC(K) 

ENDIF 

ENDIF 

ENDIF 

520  CONTINUE 
CLOSE (12) 

420  CONTINUE 

OPEN ( 12 , FILE* ' LUCILLE ' ) 

DO  700  K=  1 ,  NP 

IF(X(K) . GT. XMIN . AND .X (K) .LT.XMAX)  THEN 
IF(Y(K) .GT.YMIN.AND. Y(K) .LT.YMAX)  THEN 
IF (Z (K) .GT. ZMIN. AND. Z(K) .LT.ZMAX)  THEN 
IF (HEAD (K) .GT. 0 . )  THEN 

WRITE (12, ' (I5.4E15.5) ' )  K,X(K) ,Y(K) , Z (K) , HEAD (K) 

ENDIF 

ENDIF 

ENDIF 

ENDIF 

700  CONTINUE 
CLOSE (10) 

STOP 

END 

SUBROUTINE  REDUM(I) 

CHARACTER  DUM*2 

DO  100  J-l, I 

READ (10, ' ( A2 ) ' )  DUM 

100  CONTINUE 

RETURN 

END 
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PREFLOW  -  Program  to  extract  flow  info  from  FE  output  file 

January  1993 

$ FREEFORM 

COMMON/COOR/X (11560) ,Y(11S60) , Z (11560) , IBLK{9150.9) 

CHARACTER  FIL*12,OUTFIL (100) *5,DUM*2,TITLE*76,B*2,0*2,A*2,ZFL*1,XFL*1, YFL*1 
CHARACTER  OUT(5)  *1,  AW*76 
ZFL«= '  Z' 

XFL«'X' 

YFL«'Y' 

PRINT* ,* ENTER  INPUT  FILE  NAME' 

READ  ' (A12) '  ,FIL 
OPEN (10 .FILE-FIL) 

OPEN ( 12 . FILE- ' FLOWOUT ' ) 

CALL  REDUM ( 5 ) 

READ  (10,  '  (A76)  ' )  TITLE 
CALL  REDUM(S) 

12  FORMAT  (SIX, IS) 

READ(10, 12)  NP 
CALL  REDUM (1) 

READ (10 , 12)  NEL 
CALL  REDUM (1) 

READ (10 , 12)  MATNUM 
CALL  REDUM (1) 

READ (10, ' (SIX, Fll . 0) ' )  DATUM 
CALL  REDUM (10) 

CALL  REDUM (MATNUM) 

CALL  REDUM (10) 

DO  100  I-l.NP 

READ(10,*)  J, IB, X (I) , Y  (I) ,  Z  (I) 

IF(J.NE.I)  PRINT* ,' ERROR  --LOOP  100' 

100  CONTINUE 
CALL  REDUM (10) 

DO  110  I-l.NEL 

READ(10. ' (1015) ' )  J, (IBLK(I,K) ,K-1,9) 

IF(J.NE.I)  print*, 'ERROR- -LOOP  110' 

110  CONTINUE 
CALL  REDUM (11) 

DO  200  1-1 , NP 
READ(10, ' (16) ' )  J 

IF(J.NE.I)  PRINT*, 'ERR  -  LOOP  200' 

200  CONTINUE 
CALL  REDUM (7) 

READ (10, ' (4 OX, F12 . 0) ' )  AVFLOW 
CA1  REDUM (10) 

WRITE (12, ' (A76) ' )  TITLE 
WRITE (12,*)  AVFLOW 
DO  120  I-l.NEL 

READ(10.*)  J,VX,VY,VZ,VR, ALPHA, BETA 
IF(J.NE.I)  PRINT*, 'ERR  -LOOP120' 

CALL  CENTER  (J,XC,YC,ZC) 
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MAT=IBLK (1,9) 

WRITE (12.*)  I,  MAT, XC.YC.ZC.VX.VY.VZ.VR, ALPHA, BETA 
IF(J.EQ.O)  PRINT* ERROR  LOOP  120' 

120  CONTINUE 
CLOSE (10) 

CLOSE (12) 

STOP 

END 

SUBROUTINE  REDUM(I) 

CHARACTER  DUM*2 
DO  100  J«1,I 
READ(10, '  (A2) ' )  DUM 
100  CONTINUE 
RETURN 
END 


SUBROUTINE  CENTER { I, XC, YC, ZC) 

INTEGER  BLKNOD 

COMMON/COOR/  X(11560) ,Y(11S60) ,Z(11560) .BLKNOD (9150, 9) 

xx=o. 

YY=0. 

ZZ=0. 

DO  100  J=1 , 8 
XX=XX+X (BLKNOD (I. J) ) 

YY=YY+Y (BLKNOD ( I . J) ) 

100  ZZ-ZZ+Z (BLKNOD ( I, J) ) 

XC-XX/8 . 

YC-YY/8 . 

ZC.ZZ/8. 

RETURN 

END 
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$ FREEFORM 

PARAMETER  (NPSIZE-11S60 , NELSIZ-9150) 

DIMENSION  X  (NELSIZ)  ,  Y  (NELSIZ)  , Z (NELSIZ)  ,  VX (NELSIZ)  .VY(NELSIZ)  ,VZ(NELSIZ) 
DIMENSION  MAT (NELSIZ) 

CHARACTER  FIL*12 ,  OUTFIL (100)  *6 ,  DUM*2 ,  TITLE*  76  ,  B*2 , 0*2 ,  A* 2 .  ZFL*2  ,  XFL* 2 ,  YFL*2 

CHARACTER  OUT(5)*2 

ZFL-'ZF' 

XFL-'XF' 

YFL-'YF' 

FIL-'FLOWOUT' 

OPEN (10, FILE-FID 
READ (10, ' (A76) ' )  TITLE 
READC10, *)  AVFLOW 

10  READ(10,*, END-11)  K,MAT(K)  ,X(K)  ,Y(K)  ,Z(K)  ,VX(K)  ,VY(K)  ,VZ(K)  ,V1,V2,V3 
GO  TO  10 

11  CONTINUE 
PRINT* , ' #  EL-  '  ,  K 
NEL-K 

ZMAX-- 99999. 

XMAX-- 99999. 

YMAX--99999. 

ZMIN--ZMAX 

XMIN--XMAX 

YMIN--YMAX 

DO  300  I-l.NEL 

IF(Z(I) .GT.ZMAX)  ZMAX-Z (I) 

IF(X(I) .GT.XMAX)  XMAX-X(I) 

IF (Y(I) .GT.YMAX)  YMAX-Y(I) 

I F ( Z ( I ) .LT.ZMIN)  ZMIN-Z(I) 

IF(X(I) .LT.XMIN)  XMIN-X(I) 

IF(Y(I) .LT.YMIN)  YMIN-Y(I) 

300  CONTINUE 
PRINT*, 'X' , XMIN, XMAX 
PRINT*, 'Y' , YMIN, YMAX 
PRINT* , ' Z ' , ZMIN, ZMAX 
"  WRITE  Z- PLANES 
J- (1 . 001*ZMAX/S0) +1 
DO  400  I-1,J 
E-  (1-1) *50. +25 
IE- (E+.001) 

WRITE (OUTFIL ( I) ,'  (A2, 14) ' )  ZFL.IE 
READ (OUTFIL (I) , ' (6A1) ' )  (OUT(L) ,L-1,6) 

DO  4  L-1,6 

IF (OUT (L) . EQ . '  ')  OUT (L) »' 0' 

4  CONTINUE 

WRITE (OUTFIL (I) , ' (6A1) ' )  (OUT(L) ,L-I,6) 

PRINT*, 'WRITING  FILE  '.OUTFIL (I) 
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OPEN ( 12 , FILE=OUTFIL ( I ) ) 

WRITE (12 . ' (A76) ' )  TITLE 
WRITE (12,*)  AVFLOW 
DO  500  K=1 , NEL 

IF (Z (K) . GT.E-1 . 0 .AND. Z (K) .LT.E+1.)  THEN 
VR=SQRT(VX(K) **2+VY(K) **2+VZ (K) **2) 

WRITE (12, '  (215, 6E11.4)  ' )  K, MAT (K) , X (K) , Y (K) ,VX(K) ,VY(K) ,VZ(K) ,VR 

ENDIF 

500  CONTINUE 
CLOSE (12) 

400  CONTINUE 
"  WRITE  Y- PLANES 
J«(1.001*YMAX/50) +1 
DO  410  1=1, J 
E=(I-l)*50.+75. 

IE= (E+0 . 001) 

WRITE (OUTFIL ( I ) , ' ( A2 , 14 ) ' )  YFL, IE 
READ (OUTFIL ( I ) , ' (6A1) ' )  (OUT(L) ,L=1,6) 

DO  6  L=l, 6 

IF (OUT (L) .EQ. '  ')  OUT ( L ) = ' 0 ' 

6  CONTINUE 

WRITE (OUTFIL (I) , ' (6A1) ' )  (OUT(L) , L=l, 6) 

PRINT* ,' WRITING  FILE  '.OUTFIL (I) 

OPEN ( 12 , FILE=OUTFIL (I) ) 

WRITE (12, ' (A76) ' )  TITLE 
WRITE (12,*)  AVFLOW 
DO  510  K=1 , NEL 

IF (Y (K) . GT. E-l . 0 . AND . Y (K) .LT.E+1. )  THEN 
VR=SQRT (VX (K) **2+VY(K) **2+VZ (K) **2) 

WRITE ( 12 , ' (215 , 6E11 . 4 ) ' )  K, MAT (K) , X (K) ,Z(K) ,VX(K) ,VY(K) ,VZ(K) , VR 
ENDIF 

510  CONTINUE 
CLOSE (12) 

410  CONTINUE 

"WRITE  X- PLANES 

J= (1 . 001*XMAX/142 . 8) +1 

DO  420  1=1, J 

E=  (1-1) *142 . 8  +  178 . 5 

IE=-E+0 . 001 

WRITE (OUTFIL (I) ,' (A2, 14) ' )  XFL, IE 
READ (OUTFIL (I) , ' (6A1) ' )  (OUT(L) ,L=1,6) 

DO  8  L=1 , 6 

IF(OUT(L) .EQ. '  ')  OUT (L) = ' 0 ' 

8  CONTINUE 

WRITE (OUTFIL (I) , ' (6A1) ' )  (OUT(L) ,L=1,6) 

PRINT*, 'WRITING  FILE  ',OUTFIL(I) 

OPEN ( 12 , FILE=OUTFIL ( I ) ) 

WRITE <12, ' (A76) ' )  TITLE 
WRITE (12,*)  AVFLOW 
DO  520  K=1 , NEL 

IF (X (K) . GT . E - 1 . 0 . AND . X ( K ) .LT.E+1. )  THEN 
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VR=SQRT(VX (K) **2+VY(K) **2+VZ (K) **2  ) 

WRITE (12,  '  (215 , 6E11 . 4 )  ' )  K, MAT (K)  , Y (K)  , Z (K)  , VX (K)  , VY (K)  , VZ (K)  , VR 

ENDIF 

520  CONTINUE 
CLOSE (12) 

420  CONTINUE 

OPEN ( 12 , FILE= ' FLO ' ) 

DO  700  K=1 , NEL 

IF (Z (K)  . GT. ZMIN. AND . Z (K)  .LT.ZMAX)  THEN 
VR=SQRT(VX(K) **2+VY(K) **2+VZ(K) **2) 

WRITE (12, ' (IS , 7E11 . 4 ) '  )  K, X (K) , Y (K) , Z (K) ,VX(K) ,VY(K) ,VZ(K) , VR 

ENDIF 

700  CONTINUE 
CLOSE (10) 

STOP 

END 

SUBROUTINE  REDUM(I) 

CHARACTER  DUM*2 
DO  100  J-1,I 
READ(10, ' (A2) ' )  DUM 
100  CONTINUE 
RETURN 
END 
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FLOVEC  -  Program  to  draw  flow  velocity  vectors 
January  1993 

$ FREEFORM 
$DEBUG 

DIMENSION  X (500) , Y (500) ,V(500,3) , VMAX (3) , VMIN (3 ) ,VDIF(3)  ,  W (500, 3 ) , MAT  (500) 
DIMENSION  VEL  ( 6 )  ,  WEL  ( 6 ) 

CHARACTER  FIL* 12 , PAIR* 2 , BCK*1 , FILNAM*13 , F ( 13 ) *  1 , FILN*14 , T (6 ) *19,DOS*3 

CHARACTER  TITLE* 76 

DATA  BCK  /* W 

OPEN (3 9, FILE=' SCRATCH' ) 

5  CONTINUE 
CLINE=125. 

DO  16  1=1,6 
VEL (I) =0. 

16  WEL  (I)  =0  . 

NUM-0 
NNUM=0 
FIL* ' 

18  PRINT* ,' ENTER  FILE  NAME  OR  -  DOS  -  TO  ENTER  A  DOS  COMMAND' 

READ  ' (A12 ) ' , FIL 

IF (FIL . EQ . ' D  ' .OR.FIL.EQ. 'DOS  ')  THEN 

PAUSE  'ENTER  A  DOS  COMMAND' 

PRINT* , ' 

PRINT*, ' PRESS  RETURN  TO  CONTINUE  WITH  PROGRAM' 

CALL  UPAUSE 
GO  TO  5 
ENDIF 

READ (FIL, ' (12A1) ' )  ( F ( I ) , 1=1 , 12 ) 

F ( 13) = ' \ ' 

WRITEfFILNAM, ' (13A1) ')  (F(I) ,1=1,13) 

FILNAM=FIL//BCK 

OPEN ( 10, FILE=FIL, STATUS- 'OLD'  ,  ERR=18) 

READdO,  '  (A76)  '  )  TITLE 
READ (10,*)  AVFLOW 
1  =  1 

10  READ(10, *,END=11)  KK, MAT ( I) , X ( I > , Y (I ) , V ( I , 1) , V (I , 2) , V (I , 3) , VR 
1  =  1  +  1 

GO  TO  10 

11  N=I - 1 
CLOSE (10) 

DO  89  1=1, N 
DO  89  J=1 , 3 

89  W(I,  J)  =V(I,  J) 

IF(F(1)  .EQ. ' Z' )  PAIR* ' XY' 

IF(F(1)  .EQ. 'X' )  PAIR* ' YZ ' 

IF(F(1)  .EQ. ' Y' )  PAIR* ' XZ' 

PRINT* ,' ENTER  FACTOR  FOR  VECTOR  LENGTHS' 

READ* , FAC 

PRINT*, 'ENTER  POOL  HEIGHT' 

READ* , POOL 
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IF ( FAC . EQ . 0 . )  FAC=0. 00001 
Nl  =  2 
N2  =  3 
N3  =  l 

IF {PAIR . EQ . ' XY' )  THEN 

Nl  =  l 

N2  =  2 

N3  =  3 

ENDIF 

IF(PAIR.EQ. 'XZ' )  THEN 

Nl  =  l 

N2  =  3 

N3=2 

ENDIF 

"FIND  MAX  AND  MIN'S 
XMAX=- 9999999. 

YMAX=XMAX 
XMIN=-XMAX 
YMIN=-YMAX 
VMAX ( 1 )  =  -  9  .  E2  5 
VMAX  ( 2 ) = VMAX ( 1 ) 

VMAX ( 3 ) =VMAX ( 1 ) 

VMIN(l) =-VMAX ( 1 ) 

VMIN ( 2 ) = - VMAX ( 2 ) 

VMIN ( 3 ) = - VMAX ( 3 ) 

DO  100  1=1, N 

IF(X(I) .GT.XMAX)  XMAX=X ( I ) 

IF(Y(I) .GT.YMAX)  YMAX=Y(I) 

IF (X ( I ) . LT . XMIN)  XMIN=X ( I ) 

IF ( Y ( I ) . LT . YMIN)  YMIN=Y ( I ) 

DO  100  J=1 , 3 

I F ( V ( I , J ) . GT . VMAX ( J ) )  VMAX (J) =V(I, J) 

IF (V ( I , J) . LT . VMIN ( J) )  VMIN (J) =V ( I , J) 

100  CONTINUE 
XDIF=XMAX-XMIN 
YDIF=YMAX-YMIN 
VDIF(l) =VMAX (1) -VMIN (1) 

VDIF (2 ) =VMAX (2) -VMIN (2) 

VDIF ( 3 ) =VMAX (3) -VMIN (3) 

IF (YMAX. GT.XMAX)  XMAX= YMAX 

DIV=MAX0 (ABS (VMAX(Nl) ) , ABS (VMAX(N2) ) , ABS (VMIN (Nl) ) , ABS (VMIN(N2) ) ) 
XMAX  =  3200. 

IF(F<1)  .EQ. 'X' )  XMAX  =  2100 . 

CALL  PSTART(M, XMIN, XMAX, XMIN, XMAX) 

DO  200  1=1, N 
DO  ^00  J=l, 3 

V ( I , J) =FAC* 150 . *V ( I , J) /DIV 

200  CONTINUE 

"Z  PLANE  GRAPHICS 

IF(N1 . EQ. 1 . AND.N2 . EQ . 2)  THEN 

CALL  USET('GREE') 
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CALL  UMOVE ( 0 . , 0 . ) 

CALL  UPEN (3141.6,0. ) 

CALL  UPEN (3141. 6 ,2000.) 

CALL  UPEN ( 0 . ,2000.) 

CALL  UPEN ( 0 . , 0  .  ) 

CALL  UMOVE (856 .8, 0. ) 

CALL  UPEN (856. 8, 2000.) 

CALL  UMOVE (1356. 6, 0.  ) 

CALL  UPEN (1356. 6, 2000.) 

CALL  UMOVE (1785. ,0.) 

CALL  UPEN (1785. ,2000. ) 

CALL  UMOVE (2284 .8,0.) 

CALL  UPEN (2284. 8, 2000.) 

CALL  UMOVE (856. 8, 900. ) 

CALL  UPEN (2284. 8, 900.) 

CALL  UMOVE (856. 8, 1050.) 

CALL  UPEN (2284 .8, 1050. ) 

ENDIF 

"Y- PLANE  GRAPHICS 

IF (N1 . EQ . 1 . AND . N2 . EQ . 3 )  THEN 

CALL  USET ( ' GREE ' ) 

CALL  UMOVE ( 0 . , 0 . ) 

CALL  UPEN (3141.6,0.) 

CALL  UPEN (3141. 6, 1200. ) 

CALL  UPEN (2284. 6, 1200.) 

CALL  UPEN (1785. ,850.) 

CALL  UPEN (1356. 6, 850.) 

CALL  UPEN (856. 8, 1200.) 

CALL  UPEN ( 0 . , 1200. ) 

CALL  UPEN ( 0 . , 0. ) 

CALL  USET ('DASH') 

CALL  UMOVE ( 0 . , POOL+8 50 . ) 

CALL  UPEN (3141 . 6 , POOL+850 . ) 
CALL  UMOVE (625 ., 1200. ) 

CALL  UPEN (1356.6,650.) 

CALL  UPEN (1785. ,650.) 

CALL  UPEN (2517. ,1200.) 

CALL  USET ( ' LINE' ) 

CALL  UMOVE (1570. , POOL+850. ) 
CALL  UPEN ( 1600 ., POOL+900 . ) 
CALL  UPEN (1540. , POOL+900 . ) 
CALL  UPEN (1570. , POOL+850.) 
ENDIF 

"X- PLANE  GRAPHICS 

IF (N1 . EQ . 2 .AND . N2 . EQ . 3 )  THEN 

CALL  USET ( ' GREE ' ) 

CALL  UMOVE (0. ,0.) 

CALL  UPEN (2000. ,0.) 

CALL  UPEN (2000. , 1200. ) 

CALL  UPEN ( 0 . , 1200. ) 

CALL  UPEN ( 0 . , 0 . ) 
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CALL  UMOVE { 0 . ,850. ) 

CALL  UPEN { 500 . , 850. ) 

CALL  UPEN { 800. , 1200. ) 

CALL  UPEN { 900 ., 1200 . ) 

CALL  UPEN { 900 650 . ) 

CALL  UPEN (1050. ,650. ) 

CALL  UPEN { 1050. , 1200. ) 

CALL  UMOVE (900. ,850.) 

CALL  UPEN<1050. ,850. ) 

CALL  USET { ' DASH' ) 

CALL  UMOVE (0. .POOL+850. ) 

CALL  UPEN (2000 .  ,  POOL+850 . ) 

CALL  USET ('LINE') 

CALL  UMOVE (200. , POOL+850. ) 

CALL  UPEN  (230  .  ,  POOL+900  . ) 

CALL  UPEN (170. , POOL+900 . ) 

CALL  UPEN  (200.  ,  POOL-t-850 . ) 

ENDIF 

DO  300  1=1, N 
XI  =  I 

IF (MAT ( I ) . EQ . 3 )  GO  TO  300 
CALL  USET('BLAC') 

IF(MAT(I) .EQ. 5)  CALL  USET (' GREE ' ) 

IF (MAT (I)  . EQ . 4 )  CALL  USET  < ' MAGE ' ) 

IF (V ( I , Nl) *  *  2  + V ( I , N2 ) *  *  2 . LT .100.)  THEN 

IF ( V ( I , N1 ) . EQ . 0 . AND . V ( I , N2 ) . EQ . 0 . )  GO  TO  300 

AA=50 . * V ( I , Nl ) /MAX0 (ABS (V ( I , Nl) ) ,ABS (V(I,N2) ) ) 

BB=50 . *V ( I , N2 ) /MAX0 ( ABS (V ( I , Nl ) ) ,ABS (V(I,N2) ) ) 

V (I , Nl ) =AA 
V ( I , N2 ) =BB 
ENDIF 

TH=35 . /57 . 296 

IF(F(1)  .EQ. 'X' )  TH=4  5 . /57 . 296 
CC=COS (TH) 

SS=SIN (TH) 

I F ( MAT ( I ) . EQ . 5 )  THEN 

IF  (W  ( I ,  N3  )  .  LT  .  0  .  )  VEL  ( 5 )  =  VEL  (  5 )  +W  ( I ,  N3  ) 

IF  (W  ( I ,  N3 )  .  GT  .  0  .  )  VEL  ( 6 )  =VEL  ( 6  )  +W  ( I ,  N3  ) 

IF  (W  ( I ,  Nl )  .  LT .  0  .  )  VEL  ( 1 )  =VEL  ( 1 )  +W  ( I ,  Nl )  *SS- W  ( I ,  N2  )  *CC 
IF  (W  ( I ,  N2  )  .  LT .  0  .  )  VEL  ( 3  )  =VEL  ( 3 )  +W  ( I ,  N2  )  *SS+W  ( I ,  Nl)  *CC 
IF  (W  ( I ,  Nl )  .  GT .  0  .  )  VEL  (2 )  =VEL  ( 2 )  +W  ( I ,  Nl)  *SS- W  ( I ,  N2 )  *CC 
IF  ( W  ( I ,  N2  )  .  GT .  0  .  )  VEL  ( 4  )  =VEL  ( 4  )  +W  ( I ,  N2  )  *SS+W  ( I ,  Nl )  *CC 
NUM=NUM+1 
ENDIF 

IF (MAT (I) .NE. 5)  THEN 

IF  (W  ( I ,  N3  )  .  LT .  0  .  )  WEL  (  5 )  =WEL  ( 5 )  +W  ( I ,  N3  ) 

IF  ( W  ( I ,  N3  )  .  GT  .  0  .  )  WEL  ( 6  )  =WEL  ( 6  )  +W  ( I ,  N3  ) 

IF  ( W  ( I ,  Nl )  .  LT  .  0  .  )  WEL  ( 1 )  =WEL  ( 1 )  +W  ( I ,  Nl )  *SS  -  W  ( I ,  N2  )  *CC 
IF  ( W  ( I ,  N2  )  .  LT.  0.  )  WEL  ( 3  )  =  WEL  ( 3  )  +  W  ( I ,  N2  )  *S  S+W(I,N1)  *CC 
IF  ( W  ( I ,  Nl )  .  GT  .  0  .  )  WEL  (2)  =WEL(2)  +W(I,N1)  *SS-W(I,N2)  *CC 
IF  ( W  ( I ,  N2  )  .  GT .  0  .  )  WEL  ( 4  )  =WEL  (4  )  +W  ( I ,  N2  )  *SS+W  ( I ,  Nl )  *CC 
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nnum=nnum+i 

ENDIF 

CALL  UMOVE (X ( I ) +V  ( I , N1 ) , Y ( I ) +V ( I , K2 )  ) 

CALL  UPEN(X(I) , Y (I) ) 

IF ( V ( I , N3 )  . GE . 0 . AND . N3 . EQ . 2 )  CALL  USET ( ' BLUE '  ) 

IF(V(I,N3) .GE.O.AND.N3.EQ.2)  IC=1 

IF ( V (I <  N3 )  . LT . 0 . AND . N3 . EQ . 2 )  CALL  USET (' RED  ') 

IF(V(I,N3) -LT.0.AND.N3.EQ.2)  IC=0 

IF(V{I,N3) . GE . 0 . AND .N3 . EQ. 3 )  CALL  USET (' RED  ') 

IF ( V (I ,  N3 )  . GE . 0 . AND . N3 . EQ . 3 )  IC=0 

IF (V  (I , N3 )  . LT . 0 . AND . N3 . EQ . 3 )  CALL  USET('BLUE'  ) 

IF ( V (I ,  N3 )  . LT . 0 . AND . N3 . EQ . 3 )  IC=1 

IF ( V (I , N3 )  . GE . 0 . AND . N3 . EQ . 1 )  CALL  USET (' RED  ') 

IF ( V { I , N3 ) . GE . 0 . AND . N3 . EQ . 1 )  IC=0 

IF ( V ( I , N3 )  . LT . 0 . AND . N3 . EQ . 1 )  CALL  USET('BLUE') 

IF(V(I,N3) .LT. 0.AND.N3 .EQ. 1)  IC=1 

IF(IC.EQ.O)  THEN 

CALL  UCRCLE (X (I) ,Y(I)  ,6.) 

CALL  UCRCLE (X(I) ,Y{I)  ,9.) 

ENDIF 

IF(IC.EQ.l)  THEN 

CALL  UMOVE(X(I) -10. #Y(I) -10.  ) 

CALL  UPEN (X(I) +10. ,Y(I) +10.  ) 

CALL  UMOVE(X(I) -10. ,Y(I) +10.  ) 

CALL  UPEN (X (I) +10. , Y(I) -10.  ) 

ENDIF 

300  CONTINUE 

CALL  UMOVE (50. ,2800.) 

FILN=FILNAM 
FPDAY  =86400. 

IF(NUM.GT.O)  THEN 
IF(NNUM.GT.O)  THEN 
IF(F(1) .EQ. ' Z' )  THEN 
T (1 ) = ' AV  VEL  ( -N)  HORZ  V 
T ( 2 ) = ' AV  VEL  (+N)  HORZ  V 
T ( 3 ) = ' AV  VEL  ( -L)  VERT  V 
T ( 4 ) = ' AV  VEL  ( +L)  VERT  V 
T ( 5 ) = ' AV  VEL  (-Z)  PERP  V 
T (6 ) = ' AV  VEL  (+Z)  PERP  V 
ENDIF 

IF(F(1) .EQ. 'X' )  THEN 
T (1 ) = ' AV  VEL  ( -N)  HORZ  V 
T ( 2 ) = ' AV  VEL  ( +N)  HORZ  V 
T ( 3 ) = ' AV  VEL  ( -L)  VERT  V 
T ( 4 ) = ' AV  VEL  (+L)  VERT  V 
T ( 5) = ' AV  VEL  (-X)  PERP  V 
T ( 6 ) = ' AV  VEL  (+X)  PERP  V 
ENDIF 

IF (F (1 ) . EQ. ' Y' )  THEN 
T (1 ) = ’ AV  VEL  (-N)  HORZ  V 
T ( 2 ) = ' AV  VEL  ( +N)  HORZ  V 
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T ( 3 ) = ' AV  VEL  ( -L)  VERT  V 
T!4)='AV  vel  (+l)  vert  \* 

T ( 5 ) = ' AV  VEL  ( - Y)  PERP  V 
T ( 6 ) = ' AV  VEL  ( +Y)  PERP  V 
END  IF 

"PRINT*,'  LIMESTONE  OTHER 

XMAX-3200 . 

CALL  UWINDO ( 0 . , XMAX , 0 . , XMAX ) 

CALL  UDAREA ( 1 . , 99 . , 1 . , 99 . ) 

CALL  UFLUSH 
CALL  USET('YELL') 

CALL  UMOVE(100. ,XMAX-200.) 

CALL  UPRNT1 (FILNAM, 'TEXT' ) 

CALL  UPRNT1 ( '  LIMESTONE  OTHER \ ' , ' TEXT'  ) 

CALL  UPSET ( ' PREC' , S . ) 

DO  87  1-1,6 
CALL  UFLUSH 

YP« (XMAX-200 . ) -CLINE* (1  +  1) 

CALL  UMOVEdOO.  ,YP) 

CALL  UPRNT1 (T ( I ) , ' TEXT ' ) 

CALL  UFLUSH 

CALL  UMOVE (900. , YP ) 

I F (I . EQ . 1 . OR . I . EQ . 3 . OR . I . EQ . 5 )  CALL  UPRNT1  ( ' ' TEXT' ) 

IF ( I . EQ . 2 . OR . I . EQ . 4 . OR . I . EQ . 6 )  CALL  UPRNT1 ( '  +\ '  ,  '  TEXT'  ) 
CALL  UPRNT1  ( ABS  ( FPDAY*VEL  ( I )  /NUM)  ,  'REAL'  ) 

CALL  UFLUSH 

CALL  UMOVE (1500. ,YP) 

IF (I . EQ . 1 . OR . I . EQ . 3 . OR . I . EQ . 5 )  CALL  UPRNT1 ( ' - \ '  ,  '  TEXT'  ) 

I F ( I . EQ . 2 . OR . I . EQ . 4 . OR . I . EQ . 6 )  CALL  UPRNT1 ( '  +\ '  ,  '  TEXT'  ) 
CALL  UPRNT1  (ABS  ( FPDAY* WEL  ( I )  /NNUK)  ,  '  REAL'  ) 

CALL  UMOVEdOO.  ,100.  ) 

CALL  UFLUSH 
87  CONTINUE 
ENDIF 
ENDIF 

1222  FORMAT ( IX , A19 , F12 . 5 , 5X , F12 . 5 , 5H  FEET/DAY) 

CALL  UFLUSH 
CALL  UPAUSE 
CALL  PEND 

CLOSE (39, STATUS*' DELETE' ) 

I F ( 2 . EQ . 2 )  GO  TO  5 

STOP 

END 

SUBROUTINE  PSTART (M, XMIN, XMAX  , YMIK,  YMAX) 

COMMON/ JUNK / PCT 

CHARACTER  DEVC*4 

CHARACTER  DEV* 4 

CHARACTER  COLOR*4 , CLRBCK*4 

AA=XMIN 

AA=YMIN 

M=1 
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CLRBCK=' WHIT' 

COLOR = ' RED  ' 

PCT=0 . 45 

"  . . OUTPUT  TO  SCREEN  OR  PLOTTER . . . 

IF(M.EQ.l)  DEV= ' IBMH' 

IF (M.  EQ . 2)  DEV= ' HP4  ' 

CALL  UDEVIC (DEV) 

IF(DEV.EQ. ' IBMH' )  CALL  UDEVICE  ('EGA  ') 
IF(DEV.EQ. 'HP4  ')  CALL  UCOMPT ( 1 ) 

CALL  USTART 
CALL  UERASE 
CALL  UBACKG (CLRBCK) 

CALL  USET (' SMALL' ) 

IF(DEV.EQ. ' IBMH' )  .CALL  USET(COLOR) 

CALL  USET ( ' PERCEN ' ) 

IF(M.EQ.l)  CALL  UDAREA(1. ,99. ,1. ,99.) 
IF(M.EQ.2)  CALL  UDAREA (1 . , 75 . *PCT, 10 . , 75*PCT) 
CALL  UWINDO (0 . , 1 . 0*XMAX, 0 . , 1 . 0*YMAX) 

CALL  UOUTLN 

"CALL  USET ( ' SOFTWARE' ) 

CALL  USET ('HARD' ) 

PSFT=1 . 

IF (M. EQ . 2 )  PSFT= . 5 
CALL  UPSET (' VERT' , 100. ) 

CALL  UPSET ( ' HORI ' , 100 . ) 

IF (M. EQ . 2 )  CALL  USET ('SOFT' ) 

RETURN 

END 

SUBROUTINE  PEND 
CALL  UFLUSH 
CALL  UEND 
RETURN 
END 
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PHREAT  -  Program  to  extract  phreatic  surface  info  from  FE  file 

January  1993 

$ FREEFORM 
$ DEBUG 

COMMON/COOR/X  (11560)  ,Y(11560)  ,Z(11560)  ,XX(1000)  .YY(1000)  ,ZZ(1000) 

CHARACTER  FIL*12 , OUTFIL (100)  *5, DUM*2 ,  TITLE* 76  #  B*2 , 0*2  ,  A*2 ,  ZFL*1 ,  XFL*1 ,  YFL*1 

CHARACTER  OUT ( 5)  *1 ,  AW*76 

ZFL-'Z' 

XFL-'X' 

YFL-'Y' 

PRINT* , ' ENTER  INPUT  FILE  NAME ' 

READ  ' (A12) ' , FIL 
OPEN ( 10 , FILE-FIL) 

OPEN ( 12 , FILE- ' PHREOUT '  ) 

CALL  REDUM ( 5 ) 

READ (10, ' (A76) ' )  TITLE 
CALL  REDUM (5) 

12  FORMAT  (SIX, IS) 

READ (10 , 12)  NP 
CALL  REDUM ( 1 ) 

READ (10,12)  NEL 
CALL  REDUM (1) 

READ (10, 12)  MATNUM 
CALL  REDUM (1) 

READ (10, ' (SIX, Fll . 0) ' )  DA^UM 
CALL  REDUM (10) 

CALL  REDUM (MATNUM) 

CALL  REDUM (10) 

DO  100  1=1, NP 

READ(10,*)  J, IB, X ( I ) , Y ( I )  Z(I) 

IF(J.NE.I)  PRINT* ,' ERROR  --LOOP  100' 

100  CONTINUE 
CALL  REDUM (10+NEL) 

CALL  REDUM (11) 

DO  200  1=1, NP 

READ  (10,  '  ( 16  ,  F13 . 0 ,  F13 . 0 ,  ?.3 ,  A2  ,  A2 , 3F13 . 0)  ')  J,H,H,B,0,A,X(I),Y(I),Z(I) 
IF(J.NE.I)  PRINT*, 'ERR  -  i,OOP  200' 

200  CONTINUE 
J-l 

DO  55  1=1 , NP 

IF ( Z (I ) . NE . 0 . )  THEN 

XX (J) =X(I) 

YY ( J) =Y ( I ) 

ZZ(J) =Z(I) 

J-J+l 
END  IF 

S5  CONTINUE 
DO  77  1=1, J-l 

77  WRITE ( 12 , * )  I,XX(I) ,YY(I) ,ZZ(I) -555.0 
CLOSE (10) 
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CLOSE (12) 

STOP 

END 

SUBROUTINE  REDUM(I) 
CHARACTER  DUM»2 
DO  100  J=l, I 
READ(10, ' (A2) ' )  DUM 
100  CONTINUE 
RETURN 
END 
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